Technical Report CHL-98-24 


September 1998 


US Army Corps 

of Engineers 
Waterways Experiment 
Station 


Coastal Inlets Research Program 
Directional Irregular Wave Kinematics 


by Christopher H. Barker, Rodney J. Sobey, University of California at Berkeley 


Approved For Public Release; Distribution Is Unlimited 


am aN 


a 
wou Prepared for Headquarters, U.S. Army Corps of Engineers 


Vo. CHI : 
Under Scour at Inlet Structures Work Unit 32932 
1G = 24 


The contents of this report are not to be used for advertising, 
publication, or promotional purposes. Citation of trade names 
does not constitute an official endorsement or approval of the use 
of such commercial products. 


The findings of this report are not to be construed as an 


official Department of the Army position, unless so desig- 
nated by other authorized documents. 


B) vanres ON RECYCLED PAPER 


MBL/WHO!I 


A 


Coastal Inlets Technical Report CHL-98-24 
Research Program September 1998 


Directional Irregular Wave Kinematics 


by Christopher H. Barker, Rodney J. Sobey 


University of California at Berkeley 
Berkeley, CA 94720 


Final report 
Approved for public release; distribution is unlimited 


O 0301 0034586 4 


Prepared for U.S. Army Corps of Engineers 
Washington, DC 20314-1000 


Under Scour at Inlet Structures Work Unit 32932 
Monitored by U.S. Army Engineer Waterways Experiment Station 


3909 Halls Ferry Road 
Vicksburg, MS 39180-6199 


[ia 
US Army Corps 


of Engineers 
Waterways Experiment 
Station 


— AG fF 
8 a eon 
oy S 
fE3 ; Reet 
SS LABORAT! 


FOR INFORMATION CONTACT: 

PUBLIC AFFAIRS OFFICE 

U.S. ARMY ENGINEER 

WATERWAYS EXPERIMENT STATION 
3908 HALLS FERRY ROAD 
VICKSBURG, MISSISSIPP! 39180-6199 
PHONE: (661) 634-2502 


AREA OF RESERVATION ¢ 2.7 sqion 


Waterways Experiment Station Cataloging-in-Publication Data 


Barker, Christopher H. 

Directional irregular wave kinematics / by Christopher H. Barker, Rodney J. Sobey ; prepared for U.S. 
Army Corps of Engineers. 

183 p. : ill. ; 28 cm. — (Technical report ; CHL-98-24) 

Includes bibliographic references. 

1. Ocean waves — Measurement. 2. Wind waves — Measurement. 3. Kinematics. |. Sobey, R. J. Il. 
United States. Army. Corps of Engineers. Ill. U.S. Army Engineer Waterways Experiment Station. IV. 
Coastal and Hydraulics Laboratory (U.S. Army Engineer Waterways Experiment Station) V. Title. VI. 
Series: Technical report (U.S. Army Engineer Waterways Experiment Station) ; CHL-98-24. 

TA7 W34 no.CHL-98-24 


Contents 


1 Introduction and Review 
of Current Methods 
POTS ackpro ui Geeiysmtek sects 2 th Leb eee Ne EN? oe Nis pnt eet a 2) 22 
1.2. Methods for the Interpretation of 
terecular Waveshecords: 4/72). eons anal Ue Mine EC Roe cg Ia tae 
POG lobaleMethodsrea 2. emetic. Re ce Re 6 rem chen INGE 2 nae 
ii 2 onleocalelethiod sitet. Saar Rha) WER Pwr i ge eR ae 


2 Two Dimensional LFI Theory 
2 Problemebormulations si Seen ieerese sees meee eth Mane nte ns 2) hacer 
PPL) VAD ATA ATESY 2" She A, MERLOT eck Ch as Re A ae: 1 RAR hs td UL ce 
2.2 SE inding the Solution) yr epee. 2 5a ee Ere Mee ne earner cn 


3 LFI Method for a Point Pressure Record 

Sale hormulationsonsolutionsiee nay ae nena cea a a gen 
Sp ehormulationiomtune) Optimization eee men mee se anne 
33a iComputationy Methods stare alee? Ayick eee ae ema ne RICE oe 

d2oe ly Bre-Rrocessing oivecord’ ts | ein tat = beeen Ms Ge Pee cles 

Sore OPulnaiZationy lar OCeCUne manne eaten eee 
St el heoretieal Records: 32.88% jou bau t-te Re eee 

VAL Clroosmae laren Oi Solin o 55002000006 0000 
Stone Laboratory se Measurementsmrars ee eee te ee eee ne 
MOR IDISCUSSTOMA RAR. Te ie nS cbelc VAIL Le Ayah UR 2) a ie 


4 Three Dimensional LFI Theory 
41 eeThreedimensionaliSeast Ves 2) o.cu. pees en es Pa as 


Acs ile Backprounid asyrakt, Bybee saa Reese remen: fai an rie oS ot 
Aa5. 2 9) yar se PS ahh ere Vee ey Ee REAL Seon eee es ce 
4.4 A Local Two-Intersecting-Wave Theory................. 


iii 


A ANS Kinematics tae eas © cole Lei ce tc Seca 


5 Array of Water Surface Traces 


Heil 
ie, 
5.3 


5.4 


5.9 


5.6 


Rojsanplenareray Ch CNONNOMN 5.659 60 66 5 0,6 oa 9 wa oc 
Hormulationiof the: Optimization lessee nem ane 
Computations Methodsiu aye eu aisiey men nena 
5duly  Bre-ProcessinguomhRecords, yeast hacen ten ee 
5rot2 ee Optimizations brocedtke sien alla meee 
TheoreticalSRecordsiys oye = cps ee. Se ms eee 
Sra oimpleVavesINecordss aie) c)) aemee chen sasme nnn 
5.4.2 Records of Two Intersecting Waves ......... 
aboratony, Measurements ieaen Cara ace once 
Heo wllbaboratony, 2 xperiin en tien) ue nen ae amen 
505.2) sLaboratory, Nesultsie eee ces cases eae oe 
ASCUSSTOM OCR. (ec!) sy ee oR ana Ee eee canes ee Sree ee ee 


6 Array of Subsurface Pressure Measurements 


6.1 
6.2 
6.3 


6.4 


6.5 


6.6 


Rormulationvol Solutionmmesc: (ee oe nee waneun: | ay ey eee 
Hormillahionsofarher: Optimization men nein meine ey one 
Computation’ Methods yi -. snenone ae ee ee eee 
Godel, — WARE Pinoresetynayes Ost Invecoitels 4 5 5 oc ao go oo do & 
67322 Optimizations brocedure eatin iene ene 
heoretical Records 4. i520 Me sete) ahh eine eae cere es 
G4 Al neSingle Wave Recordswe es ws eae Snares aed: 
6.4.2 Records of Two Intersecting Waves ......... 
BielditMeasurements (ia) sca 4) 2), Seine nears Ral Dire aire bree 
6:21, GPreldi Results’ 0 24 QUE ne he te Oe Rieee eee nec 
DISCUSSION: SAWS 5) eo ss Soe ee See ee een 


7 Conclusions 


Cesk 


Pinte Workin 3 as cee sie 225) ee ren ae en Re Ameen 


Bibliography 


A Intersecting wave theory 


Al 
iNe2 
A.3 
A.4 


A.5 


Introduction’) i) 6 2 kee eg A ee ee rete 
Theoretical’ Backeround 24 54. eel lo yi een 
Analytical, Verifications y..00 pen) cues eh enone keene nen 
Numerical Verification iy) s hacen 6 ah pace ne 
VAs ee ichardsonsixtrapolationies piel sesi eee een 
Aca) aMlore'Complexiseasie)= ie) en. oe eee ee 
Depthsiof Validity 2:5 pcasapeh ath yteuey ween lan He maeime es. [ante are 


iv 


Preface 


This study was funded by the Operations and Maintenance Division, 
Headquarters, U.S. Army Corps of Engineers (HQUSACE). It is a 
product of the U.S. Army Engineer Waterways Experiment Station (WES) 
Coastal and Hydraulics Laboratory's (CHL) Coastal Inlets Research 
Program (CIRP) under the “Scour at Inlet Structures” Work Unit 32932. 
Dr. Steven A. Hughes was the CIRP Principal Investigator. 

Dr. Nicholas C. Kraus was the CIRP Technical Manager, and Mr. E. 
Clark McNair was the CHL Program Manager for CIRP. Messrs. John 
Bianco, Charles Chestnutt, and Barry Holliday were HQUSACE Program 
Monitors. General supervisory direction was provided by Mr. C. E. 
Chatham, Chief, Navigation and Harbors Division; Mr. Charles C. 
Calhoun, Jr., Assistant Director, CHL; and Dr. James R. Houston, 
Director, CHL. 


This study was conducted by Dr. Christopher H. Barker and 
Dr. Rodney J. Sobey, Department of Civil and Environmental 
Engineering, University of California, Berkeley, California, during the 
period July 1994 to September 1997. The report is a dissertation in partial 
fulfillment of the requirements for the degree of Doctor of Philosophy 
awarded to Dr. Barker by the University of California, Berkeley, 
California. 


The authors extend their gratitude to the following people for their 
assistance. Laboratory results presented in Chapter 3 were provided by 
Mr. Murray Townsend and Dr. John D. Fenton of the Australian Maritime 
Engineering Cooperative Research Center at Monash University, 
Melbourne, Australia. Laboratory data presented in Chapter 5 were 
collected with the assistance of Dr. Steven A. Hughes, Navigation and 
Harbors Division, CHL, and the technical staff at the Coastal and 
Hydraulics Laboratory, WES. Field data presented in Chapter 6 were 
provided by Mr. Gary L. Howell, Prototype Measurement and Analysis 
Branch, CHL. The FORTRAN code given in Appendix A was based on 
source code provided by Takumi Ohyama of the Institute of Technology, 


vi 


Shimizu Corporation, Tokyo, Japan. Helpful comments and review were 
provided by Professor Robert L. Wiegel and William C. Webster, 
University of California, Berkeley, California. 


At the time of publication of this report, Director of WES was 
Dr. Robert W. Whalin. The Commander was COL Robin R. Cababa, EN. 


The contents of this report are not to be used for advertising, publication, or 
Promotional purposes. Citation of trade names does not constitute an official 
endorsement or approval of the use of such commercial products. 


Chapter 1 


Introduction and Review 


of Current Methods 


1.1 Background 


Wave kinematics are involved in most processes that coastal engineers study. 
Knowledge of fluid velocities and accelerations are necessary for the study of the 
wave loading of structures through the use of the O’Brien-Morison equation. Knowl- 
edge of the kinematics near the sea bed are necessary for studies of sediment transport 
processes. For this reason, it is important that wave measurements be interpreted as 
accurately as possible. 

Despite the need for good data on the kinematics of waves, it is impractical to 
measure every parameter of interest. In a given situation, one might need to know the 
velocities, accelerations, and pressure fluctuations throughout the depth at a given 
location. While, in the laboratory, it might be possible to measure these quantities 
at many locations, it would be very expensive, and it is totally impractical in the 
field. A pragmatic approach is to measure just a few quantities, and to use a wave 
theory to compute the balance of the parameters. This approach can give a complete 
description of the kinematics from only a few measurements. Unfortunately, the 
accuracy of this approach is limited by the accuracy of the adopted theory, as well as 


the accuracy of the measurements themselves. 


A number of wave theories have been adopted in an effort to describe the kine- 
matics of waves. The most accessible and frequently used of these is Airy wave theory 
(also known as linear theory or first order Stokes theory). Airy wave theory is easy 
to use, but has a number of limitations. The source of these limitations is the sim- 
plifications of the governing equations of gravity waves that are made to linearize the 
equations, allowing a straightforward solution. These simplifications are made in the 
free surface boundary conditions and are justified by an assumption of small wave 
amplitude. Simplifying the free surface boundary conditions reduces the accuracy 
of the predicted kinematics. Unfortunately, this compromise is at the free surface, 
which is the location of the greatest fluid velocities and accelerations in waves, and 
thus frequently the most important to the forcing of structures. The assumption of 
small amplitude also renders the theory inadequate for large waves, exactly those of 
greatest interest to coastal engineers. 

In order to address these limitations, a number of high order steady wave theo- 
ries have been developed. Commonly in use are Stokes, Cnoidal, and Fourier wave 
theories. For a review of these, see Fenton (1990). In general, Stokes methods are 
successful in deep water, Cnoidal methods in shallow water, and Fourier methods in 
all depths of water. Within their limitations, all three of these methods provide ex- 
cellent predictions of the kinematics of steady waves, but are not directly applicable 
to the irregular waves commonly found in the field. With the possible exception of 
swell conditions on a very mild-sloped bottom, waves in the sea are neither steady, 
unidirectional nor monochromatic. 

Methods for the interpretation of measurements of real sea states rely on Airy wave 
theory even more heavily than do steady wave methods. With modern computer sys- 
tems, even the mest complicated of high order steady wave theories is quite accessible. 
However, such higher order theories are not directly applicable to multi-directional or 
multi-chromatic waves. The linearity of Airy theory allows superposition, in which 
any combination of waves of different frequencies or directions can be combined to 
form a solution. Superposition allows real sea states to be easily characterized by 
frequency and direction spectra. While accessible, this method results in solutions 


for the kinematics of the waves that do not satisfy the full free surface boundary 


"wave" 


Figure 1.1: Segment of record considered for global approximations 


conditions, and thus do not take into account the interactions between the individual 


waves. 


1.2 Methods for the Interpretation of 


Irregular Wave Records 


Methods used for the interpretation of irregular wave records fall into two general 
categories: global and local approximations. Global methods seek a solution that 
matches an entire measured record, or a single complete measured wave, from trough 
to following trough, or zero crossing to zero crossing (Fig. 1.1). These methods apply 
the same frequency and wave number (or set of frequencies and wave numbers) for 
all z (vertical variation) and t (time). 

Local methods, on the other hand, seek an approximation to each small local 
segment of a measured wave. In these methods, the frequency and wave number still 
apply for all z, but are allowed to vary with time, providing a separate solution in 


each small window in time (Fig. 1.2). 


Local Segment 


Figure 1.2: Segment of record considered for local approximations 


1.2.1 Global Methods 
Spectral Methods 


The most commonly used global method for the analysis of irregular waves is 
spectral analysis, coupled with superposition of linear waves. Spectral methods based 
on a single point measurement seek to define the sea state with a variance spectrum, 
E(w) which specifies the contribution to the variance of the water surface at each 
frequency, w. In practice, the spectrum is defined in discrete form, with the water 


surface described by the superposition of many linear waves, 


M 
aot) = Se Am COS (Km (x CoS Am + ysin On) — Wmt + Am) (ali) 


m=1 
where 77 is the elevation of the water surface, w,, and k,, are the frequency and wave 
number of the mth wave, 0,, is the direction of propagation of the mth wave and a,, 
is the phase of the mth wave at the origin of the coordinate system. 


The amplitudes are found from the discrete variance spectrum, 
Om = V 2E (wm )Aw (LZ) 


where Aw is the frequency spacing of the discrete spectrum. The frequency and wave 


number of each of the M waves are related by the linear dispersion relation, 


w? = gkm tanh kph (1.3) 


where g is the acceleration of gravity, and h is the mean water depth. The dispersion 
relation defines the phase speed of each component, allowing each separate component 
to move independently, without interaction with one another. Any Eulerian current 
is not included in this relation, and is generally ignored in spectral methods. 

The discrete variance spectrum, EF’, can be computed from a water surface or sub- 
surface pressure record through the use of variations of the Fast Fourier Transform 
(Dean and Dalrymple 1991; Newland 1993). A point measurement provides no in- 
formation about the direction of propagation of the waves, 6,,. It is usually assumed 
that all the waves are propagating in the same direction. Once the amplitude, fre- 
quency, wave number and phase of each individual wave are defined, the kinematics 
and dynamics of the wave field can be computed by superimposing the kinematics of 


each individual wave, as predicted by linear wave theory. 


Directional Spectra When measurements are taken by an array of instruments, 
the directional nature of the sea can be described by a directional spectrum, S(w, @). 
In this case, the water surface is represented by a large number of linear waves of 


different frequencies and directions: 


M N 
(2508) = es Ss Amn COS (km(x cos, + y sin O,) — Wmt + Amn) (1.4) 


m=1 n=1 


where: 
Onn N29) Com a NONG (1.5) 


and Aw and A@ are the sample spacing of the discrete spectrum in frequency and 
direction space. w,, and k,, are once again related by the linear dispersion relation 
(Eq. 1.3). The directional spectrum is usually broken down into two parts, the one 
dimensional variance spectrum, E(w), and the directional spreading function (DSF), 
D(w, 6): 

S(w, 0) = E(w)D(w, 8) (1.6) 


where E(w) is the variance spectrum defined above, and 


‘3 iP D(w,0)d0dw = 1 (1.7) 


There is a great deal of literature about how to best determine the DSF from a variety 
of arrangements of measurements. Current practice is reviewed in Mansard (1997). 
Unfortunately, all of these methods are limited in their ability to define accurately 
the DSF. 

The time series measurements used to determine the variance spectrum commonly 
include in excess of 1000 points in time. This much data in the time domain allows 
a very high resolution computation of the spectrum of a stationary process in the 
frequency domain. In the spatial domain, by contrast, there are only as many data 
points as there are instruments in the particular measurement array. This may be 
as few as three, and perhaps as many as a dozen, but it is too expensive to use 
many more than that. As a result, there is little information to define the DSF. For 
example, when data from an array of three instruments is analyzed by the standard 
linear method, the DSF can be specified by only five independent coefficients (Dean 
and Dalrymple 1991). While these few coefficients may serve well for determining 
integral properties, such as the mean direction and radiation stress, it is not enough 
information to accurately specify the complete kinematics. 

Another difficulty arises when determining the spectrum from measurements other 
than wave staffs. Most often, subsurface pressure gauges or a combination of pressure 
gauges and orthogonal velocity gauges are used. In this case, the measured quantity 
must be related by a transfer function to the equivalent water surface. The transfer 
function is most commonly determined from linear wave theory. This use of linear 
wave theory once again contributes to errors in situations where linear theory is not 
entirely appropriate. Bishop and Donelan (1987) discuss the difficulties in determin- 
ing the transfer function from a subsurface pressure measurement to the equivalent 
water surface. 

The commonly used methods for determining the DSF are statistical methods, 
in that they rely on the assumption that the phases of the individual components 
are randomly distributed. This assumption allows the DSF to be computed by dis- 
regarding the phase information. Unfortunately, without the phase information, it is 
impossible to reconstruct the detailed kinematics, only statistical descriptions can be 


formulated. 


Even if a frequency or frequency-direction spectrum has been accurately deter- 
mined, these methods have a number of shortcomings when used to predict the kine- 
matics of irregular waves. Shortcomings include the inaccuracies inherent in linear 
wave theory, particularly in shallow water and with large waves. 

An additional difficulty arises from the superposition of many waves. This problem 
is known as high frequency contamination of the crest kinematics (Forristall 1985; Lo — 
and Dean 1986; Bishop and Donelan 1987). Fundamentally, the difficulties arise 
from the approximations made by linear wave theory in the free surface boundary 
conditions. In linear theory the free surface boundary conditions are applied at the 
mean water level, and thus predictions made above that level are strictly out of 
the solution domain. If the full free surface boundary conditions are not satisfied, 
the resulting predictions will be inaccurate, particularly near the free surface. In 
particular, the hyperbolic function quotients that define the vertical variation of the 
kinematics become very large in the region above the MWL for the high frequency 
(and high wave number) components. This results in substantial high frequency 
fluctuations in the predicted kinematics near the crest. 

In an attempt to reduce the inaccuracies in linear superposition’s predictions of the 
near surface kinematics, empirical modifications to Airy theory have been adopted 
(Wheeler 1969; Lo and Dean 1986). This method, known as Wheeler stretching, 
locally adjusts the vertical dimension to prevent the evaluation of the hyperbolic 
quotients from being evaluated above the MWL. The result is a hybrid global-local 
method in that the frequency, wave numbers, and amplitude of each wave are deter- 
mined globally, but the coordinate system is defined locally, varying with time. 

The stretching method produces predictions of crest kinematics that seem to 
match measured data better than simple linear superposition, but it no longer satis- 
fies either the Laplace equation (mass conservation), or the full free surface boundary 
conditions. 

Another attempt to improve on the accuracy of determining the kinematics from 
directional spectra is the recent work by Prislin et al. (1997), and Prislin and Zhang 
(1997). This method seeks to reconstruct wave kinematics from a directional spec- 


trum using a second order Stokes-type interacting wave theory. The spectrum is 


decomposed into a set of individual free waves, with from one to five directional free 
waves per frequency. The effects of the second order interactions are subtracted from 
the measured record to determine the amplitudes and phases of the individual free 
waves through an iterative procedure. This results in an expression for the potential 
function and water surface that includes a full set of many free waves, and the corre- 
sponding second order bound modes. The method succeeds in globally reproducing 
the measured kinematics in the given deep water field records quite well, although 
the largest errors are in the vicinity of the crest, where velocities are greatest, and 
accuracy is most important. 

Another limitation of the Prislin and Zhang method is that the nonlinear interac- 
tions are computed through the use of a Stokes-type perturbation expansion in wave 
steepness. Based on experience with steady waves, a second order expansion of this 
type is likely to be adequate in deep water, but if the method were to be applied 
in transitional, and especially in shallow water, a much higher order representation 
would be necessary (Fenton 1990). While it is theoretically possible to extend the 
method in this manner, it would increase the magnitude of the computation substan- 
tially, perhaps prohibitively. Representing an entire record in the global sense requires 
many interacting free waves. Considering their interactions at high order would pro- 
duce huge numbers of interacting terms that must be considered. An alternative, 
local, approach would need to consider far fewer interacting waves. 

Other global methods rely on zero crossing analysis to identify particular waves 
that are then analyzed by using steady wave theory for a wave of the same height and 
period. This approach can provide an order of magnitude estimate for the kinematics, 
but does not take into account the detail of the record, and thus can not be expected 
to consistently provide better than order of magnitude accuracy. 

To include the detail of the record in wave by wave analysis, Dean (1965) adapted 
his numerical stream function method to irregular waves, seeking a Fourier expansion 
for the stream function that includes both cos j(kx — wt) and sin j(kx — wt) terms. 
The method optimizes the Fourier amplitudes to best solve the free surface boundary 
conditions at the measured water surface elevations from trough to following trough. 


While this approach takes into account the detail of the record, the method includes 


z = (water surface) 


Figure 1.3: Coordinate system for Lambrakos-Baldock-Swan method. 


only a single free mode, with all other included components being bound modes 
traveling with the wave at the same phase speed, thus representing an asymmetric 
wave of permanent form. A solution that does not allow a change in form is unlikely 
to accurately capture waves in deep water, where frequency dispersion can lead to 
transient extreme waves, or, indeed, in shallow water where shoaling effects cause a 
change in form as the waves progress. 

Seeking to improve on global methods, Lambrakos (1981) developed a method for 
determining the kinematics of two dimensional irregular waves that includes many 
free modes, and thus unsteady motion. Baldock and Swan later refined Lambrakos’ 
method, applying it specifically to large transient waves, first in deep water (Baldock 
and Swan 1994), and then in shallow water (Baldock and Swan 1996). Baldock 
and Swan’s method adopts the following potential function that is a double Fourier 


expansion in space and time: 


N M 
NG Ft) = yy ee cosh (nkz) (Anim cos (nkzx — mwt) + Brim sin (nkz — muwt)) 
n=1 m=1 


(1.8) 


where k, w, Anm, Brym are global constants, and z = 0 at the bed (Fig. 1.3). 
In this potential function each frequency component (mw) has a corresponding set 


of wavelengths (nk), allowing each component to travel at distinct phase speeds. 
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This approach can accommodate both steady and unsteady wave forms by including 
both a set of free waves and a corresponding set of bound modes. Periodicity over 
time and space is assumed. The record is not likely to be periodic, so the duration 
of the record is assigned as the fundamental period, (27/w), and the fundamental 
wavelength, (27/k) is found as part of the solution. The method uses a nonlinear 
optimization to find the coefficients, A,,., and By, m, the fundamental wavenumber, 
k, and the Bernoulli constant that produce a solution that globally satisfies the full 
free surface boundary conditions. 

In order to accommodate the unsteadiness of the wave profile, the evolution of 
the wave in space must be known. This is accomplished by applying the free surface 
boundary conditions at many locations in time and space. Since the elevation of the 
water surface usually is measured only at a single spatial location, it is predicted at 
the other locations as part of the solution. The potential function (Eq. 1.8) exactly 
satisfies the bottom boundary condition and mass conservation. In order to find the 
solution that best fits the free surface boundary conditions, the method seeks a set 
of coefficients that minimizes the sum of squares error in both of the free surface 
boundary conditions over a grid of nodes in time and space. 

Baldock and Swan identified a difficulty in this basic method, as the squared error 
was equally considered at all of the nodes, most of which were at locations in which 
the water surface elevation was also unknown. This resulted in solutions in which the 
error in the boundary conditions was greatest at the location of the measurement. 
The difficulty was mitigated by a weighting function, multiplying the dynamic free 
surface boundary condition errors at the measured location by a factor of 50 in the 
sum of squares calculation, to force the solution to match well at that point. 

Comparisons of their results with laboratory data were quite good, but the method 
has a number of limitations. These include a huge matrix of unknown coefficients that 
must be found simultaneously (2(/7N +1) unknowns, with typical values of M and N 
of 18, for a total of 650 unknowns), and the method must solve for the water surface far 
from the measurement location, removing the focus from the actual measured data. 
The method makes no assumptions about the steadiness of the wave field; however, it 


extends the horizontal bottom assumption and introduces unnecessary complication 
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in its attempt to capture an entire segment of an unsteady wave field with a single 
expression. This complication would become at least an order of magnitude greater 


if the method were extended to include the second horizontal dimension. 


1.2.2 Local Methods 


Nielsen (1986, 1989) introduced a method that uses a local frequency and linear 
wave theory to find the location of the water surface from a pressure record. In this 
method, the wave in the region of each measured point is considered to be a small 


segment of a linear wave, so that the pressure record is locally represented by a sine 


function. 
p = An sin (Wptn — Qn) (1.9) 
so that: 
1 02p 
n= 4{-—-— 1.10 
Bs p Ot? Gey) 


Estimating the curvature by finite differences, the local frequency becomes: 


—Pn-1 ate 2Dn — Pn+l1 
Sap pi Se ae ial 
= pAb? ay) 


where At is the time spacing between points. Alternatively, w can be computed from 
trigonometric identities, 

On = 500s" (Batre | (1.12) 
which is exact for a sinusoid. 

Once the frequency has been determined, the wave number is computed from the 
linear dispersion relation, Eq. 1.3. The water surface elevation can then be computed 
with the linear pressure response function: 

_ Pn cosh(knh) 


= 113 
pg cosh(ky Zp) Was) 


Nn 


where z, is the elevation of the pressure measurement above the bed, and h is the 


mean water depth (As the method is a local method for the interpretation of a point 


12 
measurement, the local mean water depth is used, rather than the global still water 


depth). This may be adapted somewhat by using a variation of Wheeler stretching: 


cosh (k,, {h + 2 
Tn = poe) (1.14) 


Nielsen also presented an alternative approach, which uses a semi-empirical trans- 
fer function derived from Fourier wave theory to compute the water surface elevation 
from the measured pressure and local curvature of the pressure record. While quite 
accurate at reproducing the water surface, neither method supplies the kinematics, 
nor do they satisfy either mass conservation and or the free surface boundary con- 
ditions. Despite these limitations, the efficacy of these methods demonstrates the 
potential for the local approach. 

To interpret bottom pressure measurements in the context of the kinematics while 
preserving the full governing equations, Fenton (1986) presented a method that em- 
ploys a local polynomial approximation to the complex potential function. In Fenton’s 
method, the potential function and water surface are represented by separate poly- 


nomials in each small window in time: 


M 
o(x — ct,y) + iv(ae — ct,y) = iene a (1.15) 
M 
Nee) = SPE — ct) (1.16) 
j=0 


where z = x+y, y = 0 at the bed, 77 is the water surface, c is the wave celerity, 
and the a; and 6; are real. The wave is assumed to propagate at speed, c, without 
change in form. While steadiness is not a valid assumption in the global sense, it is 
only applied locally, within a small window in time. 

The method solves for the coefficients, a; and b;, and the wave celerity, c, that 
satisfy the full nonlinear free surface boundary conditions and fit the measured pres- 
sure record. This approach provides the complete kinematics and satisfies the full 


governing equations. Based on a polynomial variation with depth, it works well in 
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shallow water, where Cnoidal methods with polynomial vertical variation are theo- 
retically appropriate. At the same time, it may not be applicable in transitional or 
deep water, where Stokes type methods, based on an exponential variation in the 
vertical, are theoretically more appropriate. In a later paper, Fenton and Christian 
(1989) presented a simplified version of the method that required less algebraic ma- 
nipulation, and was still found to be effective in shallow water. In this case, the local 


wave celerity was assumed to be that given by long wave theory: 


C= SG (L107) 


While this is a reasonable approximation for long waves in shallow water, it was 
not appropriate for shorter waves in deeper water, as might be expected from the 
polynomial form and the long wave celerity. 

To find a method that could work in any depth of water, Sobey (1992) developed 
the Local Fourier Method for Irregular waves (LFI). This approach employs a po- 
tential function represented by a low order Fourier expansion in a small window in 
time. It is a method derived for the analysis of a point water surface trace. Local 
frequency, wave number, and the Fourier coefficients are sought that fit the measured 
record and the full free surface boundary conditions. The LFI method provides the 
' complete kinematics, satisfies the full governing equations, and is successful in all 
depths of water. A more complete description of this method follows in Chapter 2. 
Sobey’s method shows a great deal of promise, but was only applied to the analysis 
of a water surface measurement at a single location. 

Pressure sensors are easy to deploy, and thus are often used to measure waves, 
particularly in shallow and transitional depth water. The LFI method is extended 
in this dissertation to the interpretation of wave measurements from a point pressure 
sensor in Chapter 3. 

Point measurements provide no information about the directionality of the mea- 
sured wave field. Wind seas are not uni-directional, and knowledge of directionality 
has been shown to be very important in the prediction of the kinematics and forcing 
of structures by real seas (Dean 1977; Forristall et al. 1978). 


Arrays of measurements are frequently used in order to capture the directional 
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nature of the wave field. The LFI method has been further extended in this disserta- 
tion to the interpretation of arrays of water surface measurements in Chapter 5, and 


arrays of pressure sensors in Chapter 6. 
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Chapter 2 
Two Dimensional LFI Theory 


It is common to deploy a single measurement device to gather information about 
the wave field at a given location. Measurement at a single point in space does not 
give any information about the directionality of the sea state, but if the wave field 
is assumed to be locally two dimensional, a reasonable description of the kinematics 
can be established. This chapter outlines the LFI method as it is applied to the 
interpretation of a time series from measurements taken at a single location in space, 


such as recorded by a pressure gauge or wave staff. 


2.1 Problem Formulation 


The problem formulation for two dimensional irregular waves has much in common 
with classical steady wave theory. The flow is assumed to be incompressible and 
irrotational. The kinematics can therefore be represented by a potential function, 
o(x,z,t), where 
_ 0 

Ox 


_ 96 


U 


and u and w are the horizontal and vertical velocities, respectively. 
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Figure 2.1: Coordinate system used for the two dimensional method 


Field Equation The field equation is mass conservation for irrotational flow, in 


the form of the Laplace equation: 


Od 0d 
Wie = ey 
¢ Ox e Oz? : 2) 
Boundary Conditions The boundary conditions are the bottom boundary condi- 
tion for a locally horizontal bed (BBC), 
eee 


= apie 0 at a) (2.3) 


(20) 


the kinematic free surface boundary condition (KFSBC), 


w—-——-—u— =0 at z=N (2.4) 


1 rt 
—+=-uv+-w?+9n-B= at z= (2.5) 


where 77 is the elevation of the free surface, g is the acceleration of gravity, and B is 
the Bernoulli constant. 
The kinematic free surface boundary condition (Eq. 2.4) includes both the time 


and spatial gradients of the water surface. When working with a subsurface gauge, 


17 


the location of the water surface is unknown and part of the solution. When working 
with data from a water surface probe, the location of the water surface at discrete 
points in time is known, but the gradients in time and space are not. In order 
to accommodate this lack of data, the kinematic free surface boundary condition 
is transformed to eliminate the gradients of the water surface. Following Longuet- 
Higgins (1962), the kinematic condition is subtracted from the the dynamic condition 


differentiated following the motion: 
MKFSBC = —g(KFSBC) + = (DFSBO) (2.6) 


Resulting in a modified kinematic free surface boundary condition: 


qP Uap Uo = = (0) at 227) 


This new form of the boundary condition does not include the gradients of the water 
surface. Applying both the modified kinematic free surface boundary condition and 


the dynamic free surface boundary condition completes the formulation. 


Observational Equations In steady wave theory, periodic boundary conditions 
are also imposed, forcing the solution to be periodic in both space and time. For 
irregular waves, however, the periodicity is not known. Rather, the local solution 
is defined by a local segment of a measured record within a small window in time, 
together with the field equation and the bottom and free surface boundary conditions. 

In order to define a solution that fits the measured record, observational equations 
are identified. These equations will be different depending on which quantity has 
been measured. In the case of a water surface measurement, they are the free surface 
boundary conditions, applied at the measured elevations at a number of points in 
time throughout the window (Sobey 1992). In the case of a pressure measurement, 
they are the Bernoulli equation, applied at the elevation of the measurement, and 


also at a number of points in time within the window (see Chapter 3). 
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Solution in Each Local Window A form for the potential function in each local 
window is based on Fourier steady wave theory. This is the potential function used 
by Sobey (1992): 


@(z,2,t) =Ux + yA jE sin i(k — wt) (am) 


where U and hf are the known cera Eulerian current and mean water depth 
(see section 3.3.1 for a discussion of these important parameters), J is the truncation 
order of the Fourier series, A; are the Fourier coefficients, and w and k are the 
local fundamental frequency and wave number. The above potential function exactly 
satisfies mass conservation and the BBC. This form for the potential function is 
periodic in space and time, however the periodicities are not defined apriori, but 
found to fit the record, defining a local frequency and wave number. 

The measured record is broken down into individual segments, each in a separate 
window in time. In each window, a different set of the parameters w, k, kx, and A; 
are found to fit the segment of the measured record. This represents that segment as 
a piece of a larger, periodic wave. The entire record is then represented by separate 


potential functions, each applied to a particular window in time. 


2.1.1 Dynamics 


The potential function provides the complete kinematics, and the dynamics are 


found through the unsteady Bernoulli equation: 


Ft 5 +u?) +a24+2 -B=0 (2.8) 
where p is the mass density of the water, and p is total pressure. Total pressure below 
the water surface can be broken down into three components. Atmospheric pressure 
(pa) is the pressure of the atmosphere at the water surface, hydrostatic pressure (p,) 
is the pressure in the water column due to gravity, in the absence of motion. Dynamic 


pressure (pq) is the component due to the motion of the fluid. Separating these three 


components focuses attention on the wave motion. 


P= Pa + Pa + Ph where Ph = —pgz jor Zhe (2.9) 
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When atmospheric pressure is defined as zero, and dynamic pressure is substituted 


into Eq. 2.8, the result is: 


xo) 1 2 2 Pd — 
= = ) — — == » 
ag oC ee at (2.10) 


The Bernoulli Constant 


In a potential flow, the Bernoulli constant is the same throughout space and 
time. For the special case of wave motion, the value of the Bernoulli constant can 
be computed if the kinematics are known (Longuet-Higgins 1975). The Bernoulli 
equation (2.8) and the bottom boundary condition are applied at the bottom: 

set stot Bao (2.111) 
where wu,, pp, and z, are the velocity, pressure, and elevation at the sea bed. If the 


flow is periodic, a time average over a period results in: 
T= De a 
pu toa t > -B=0 (BI) 


where the over-bars indicate time averaging. 

In the case of steady periodic wave motion, the total vertical momentum at any 
horizontal location must be the same at the beginning and end of a period. In order 
for this to be the case, the vertical momentum averaged over a period must be a 
constant. In order for the momentum to remain constant, the time averaged net 
vertical force on the water column at that point must be zero. If the pressure at the 
water surface (p,) is taken to be zero, then the force of the pressure at the bed must 
be equal to the gravitational force on the column of water, so that the mean pressure 


on the bed is hydrostatic, 
Ps = pg(7 — 2s) (2.13) 
resulting in a simple and exact expression for the Bernoulli constant. 


ast i 
B= gi + su; (2.14) 
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With z defined to be zero at the mean water level (7) and Eq. 2.7 as the potential 
function, B becomes (Sobey 1992): 


2 
=i mele hae gkA; 5 
i Se su a7 > (5) (2.15) 


and all the terms in Eq. 2.10 are defined by the potential function, allowing the 


dynamic pressure to be computed from the kinematics. 


2.2 Finding the Solution 


The bulk of the LFI method is the process of determining appropriate values for 
each of the parameters, w, k, kx, and A,, in the potential function (Eq. 2.7) for each 
window in time. The result is a set of potential functions that completely describe the 
kinematics of the wave field in the region of the measurement, each separate window 
in time being described by a unique and independent potential function. 

In the case of subsurface measurements, the location of the water surface is also 
required. Determination of the parameters of the potential function, as well as the 
location of the water surface, is accomplished through nonlinear optimization that 
matches the measured record and minimizes the error in both free surface boundary 
conditions. Sobey (1992) applied the LFI method to the interpretation of a single 
water surface trace. The following chapter applies the method to a subsurface pressure 


record. 


Chapter 3 


LFI Method for a Point Pressure 
Record 


Pressure gauges are relatively inexpensive and easy to deploy, and as a result are 
frequently used in the field to measure waves, particularly in shallow water. The 
simplest instrument consists of a single subsurface pressure gauge. While a single 
point measurement provides no directional information, a reasonable estimate for the 
kinematics of the measured waves can be obtained. This chapter describes a technique 
_ for applying the two dimensional LFI method presented in the previous chapter to 


the analysis of a point subsurface pressure trace. 


3.1 Formulation of Solution 


The flow is taken to be two dimensional, with the governing equations described 


in Chapter 2. These include the potential function, 


hjk(h +z 
WG, 21) = er ee 2 sin j(kx — wt) (3.1) 


the modified kinematic free surface boundary condition (f*), 


22 


Od Ow 
Fg oad KD Sic) (eee aes 
jf = Fp ah ED A Pa) 
an oe AL Ow 
dr Oz 
Ou Ow 
=F a + w= =0 at B= i 
(3.2) 
the dynamic free surface boundary condition (f”), 
Oo = 
D 2 2 me puts 
i = Fe se ath Gj Ba) 2= 7 (3.3) 
and the unsteady Bernoulli equation ( f8). 
0d 1 (de 
BOON ea o BI ipa Hels rams 
if Sar al ae B= 0 (3.4) 
with the Bernoulli constant defined as: 
fh, 1 1 AlnAle \* 
B = aCe oe 33, 
a “ 4 X (5) ce) 


The unknown parameter, x, appears in the potential function only when coupled 
with the parameter, k. It is convenient to solve for the non-dimensional parameter, 
kx, essentially a phase parameter in the potential function. 

The action of waves is greatest near the surface, so that the fluid velocities and 
accelerations are largest near the surface and decay rapidly with depth, particularly 
in deep water. In addition, for pressure sensors placed at or near the bottom, the 
vertical component of the velocity is damped out completely, as expressed by the 
bottom boundary condition. The problem, then, is to extrapolate the motion of the 
fluid up to the surface, where the motion is greatest, using only data extracted from 


the subsurface, where the motion is much less. 
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Mathematically, this is an ill posed problem, as the flow is governed by an elliptic 
equation (Laplace, Eq. 2.2), in which the solution is determined by the boundaries, 
but there are only data at a single location. Some of these difficulties can be sur- 
mounted through the application of the free surface boundary conditions. While there 
are no data measured near the surface, the free surface boundary conditions remain 
appropriate, and necessary to define the solution at that boundary. 

When the LFI method was applied to a water surface trace (Sobey 1992), the 
location of the the water surface was known, and the boundary conditions could 
be directly applied at that location. Unfortunately, when working from a subsurface 
pressure record, the location of the water surface is not known. While the free surface 
boundary conditions are well defined, the fact that the location at which they must 
be applied is not known makes the problem more difficult. This is the complication 
that leads to the difficulties in finding full nonlinear solutions to all free surface flow 
problems. In order to apply the free surface boundary conditions when working with 
a subsurface record, the location of the water surface must be found, together with 
the potential function, in each window. 

In order to locate the free surface, the water surface is defined at N surface 
nodes equally spaced in time throughout the window. The elevation of these nodes 
“is unknown, and will be sought as part of the solution. Including the water surface 
nodes as part of the sought solution introduces N additional unknown parameters for 
a total of 3+ J+ N unknowns in each window (k, kz, w, Ay... Ay, ---7n). The 


free surface boundary conditions, Eq. 3.3 and 3.2, are applied at each surface node. 


3.2 Formulation of the Optimization 


Finding the unknown parameters in a nonlinear system of algebraic equations is 
known as nonlinear optimization. The system in this case consists of the nonlinear 
Bernoulli equation and the nonlinear free surface boundary conditions. The free 
surface boundary conditions are nonlinear in two respects, involving second order 
terms in dynamic free surface boundary condition (the u” terms and B in Eq. 3.3), 


and second and third order terms in the modified kinematic boundary condition (the 
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Figure 3.1: Schematic of system of equations in a window 


uz and u?2 terms in Eq. 3.2). There is also a significant nonlinearity introduced b 
Ot Or q g y vy 


the application of the boundary conditions at the unknown and varying free surface. 


Observational Equations The given form for the potential function could repre- 
sent any periodic flow, subject to the bottom boundary condition. The FSBCs define 
the flow as a gravity-constrained, free surface flow. The observational equations are 
the equations in the system that force the solution to fit the given record. For a 
subsurface pressure record, this is the Bernoulli equation, applied at the location of 
the pressure measurement. The required number of independent equations are es- 
tablished by applying the Bernoulli equation at a number of times throughout the 
window considered. The error in the Bernoulli equation is the difference between the 
measured dynamic pressure and that computed from the kinematics defined by the 
potential function. The solution is the set of parameters in the potential function 
and the set of water surface nodes that produces a predicted dynamic pressure that 
matches the measured record, while simultaneously satisfying the FSBCs. 

A system of equations is specified if there are as many independent equations as 
unknown parameters in the system. If there are more equations than unknowns, the 
solution can be defined as that which results in the smallest squared errors in the 


equations. This least squares formulation is also appropriate for a uniquely defined 
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system, with an expected error of zero in all the equations. In order to specify the 
solution to the LF] formulation, the boundary conditions (f? and f*) are applied at 
each of the water surface nodes, yielding 2V equations, and the Bernoulli equation 
(f®) is applied at I nodes on the pressure record within the local window (Fig. 3.1). 
There are a total of 3+ J+ N unknowns (k, kz, w, Ai...Ajz, and m...n). The 
system is uniquely specified for / + N =3+ J and overspecified for /+ N >3+J, 


resulting in the following least squares optimization: 


minimize O(X )= SAP 2G has ten Pat CS Tins t tn)? 
=P So F(X: lie Bop ta) (3.6) 
7=1 


XS = (6, the A ooo Alig ih ooo iy) 


where t, are the locations in time where the water surface nodes, 7,, are sought. 
t; are the times within the window where the Bernoulli equation is applied, Py; 
are the measured dynamic pressures at t;, and z, is the elevation of the pressure 
gauge. Overspecification, particularly with additional nodes on the pressure record, 
is advantageous for an actual record as a way to minimize the effect of any noise in 
the measurements. Additional nodes on the pressure record also serve to emphasize 
the measured data, which directly define the local kinematics. The details of the 


solution to this system of equations is given in the following section. 


3.3 Computation Methods 


The LFI method can be broken down into the following sequence of steps: 
1. Pre-processing of record. 


(a) Determine estimate for level of noise in the record. 
(b) Determine MWL and subtract hydrostatic pressure from record. 


(c) Determine estimate for magnitude and direction of Eulerian current. 
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(d) Define a continuous record from the discrete observations by cubic spline 


interpolation. 
(e) Specify spacing of output locations. 
(f) Compute mean zero crossing frequency. 


(g) Non-dimensionalize record and all parameters. 
2. Primary values of numerical solution parameters are chosen. 


(a) Window width (70) 
(b) Order of solution (J) 
(c) Number of sample locations on record within each window (J) 
(d) Number of water surface nodes within each window (JN) 
3. For each selected output location, a window in the record is defined, and an 
LFI solution is computed. 
(a) Initial guess for the optimization is computed. 


(b) Full nonlinear optimization for all unknown components of the potential 
function is computed. 


(c) Results are checked for spurious solution. 


i. If no solution, or a spurious solution, is found the solution parameters 


are adjusted, and the optimization repeated. 


ii. If a good solution is found, progress to the next window. 


3.3.1 Pre-Processing of Record 
Accommodating Measurement Error 


Measurement error can be a major source of difficulty, as the method relies on 
the details of a small segment of the record. In Sobey’s work (1992), the primitive 
kinematic free surface boundary condition (Eq. 2.4) was applied at the measured water 


surface, requiring an estimate for the time gradient of the surface. The estimation 
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of the gradients from a measured record is very sensitive to error in the record, and 
thus Sobey found it necessary to smooth field and laboratory measurements with a 
simple moving average filter. 

There are no pressure gradient terms in the Bernoulli equation, so the application 
of the LFI method to a measured pressure trace should be less sensitive to noise. 
Nevertheless, measurement error remains a problem in applying the LFI method toa 
pressure record. When there is obvious noise in the record, an alternative to smooth- 
ing is to substantially overspecify the system of equations. By taking many samples 
on the pressure record in each window, the least squares optimization may accommo- 
date the errors in the measurements. This is preferable to smoothing, as no smoothing 
algorithm can reproduce lost data, but will rather impose some apriori assumption 
about the nature of the record on the results. Accommodating the measurement error 
with many samples on the pressure record allows the least squares optimization to 
find the best fit without imposing any further assumptions on the results. This was 
the approach taken with the presented laboratory data (Section 3.5). Records gener- 
ated analytically by Fourier wave theory contained no error, and therefore required 


no special accommodation. 


- The Mean Water Level 


The mean water depth, h, is included in the potential function (Eq. 3.1). This 
is an unknown value when the only data available is a pressure record. While an 
approximate depth of water is likely to be known at a given deployment site, the 
water level fluctuates due to astronomical and storm tides. These processes happen 
over a much larger time scale than the periods of wind generated waves. An estimate 
for the local mean water level (MWL) can be determined by averaging the measured 
pressure (after subtracting out the atmospheric pressure) over a period of time that 
must be larger than a typical wave period, but smaller than the period of tidal 
fluctuations, to accommodate changes in the mean water level. This mean pressure 
is the local hydrostatic pressure, which is used to compute the mean water level, and 


is subtracted from the pressure record to define the dynamic pressure. While this is 
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not an exact procedure, the MWL is not a critical parameter, as h is only used to 
locate the origin of coordinates in the potential function. Note that the result is a 
local MWL, which could be quite different than the global still water level (SWL). 
The mean water level is the more appropriate reference frame, as the method is local, 
and only considers a small window in time at once. The elevations of the water 
surface nodes are found independently, without any requirement that the resulting 
mean water level be at the origin. Thus it is not necessary that the origin is at the 
exact MWL, only that the it is near the surface. In shallow water, the bed could be 
used as the origin (Fenton 1986), but this would not work well in deeper water, as 


the area of primary interest, the water surface, would be far from the origin. 


The Eulerian Current 


The depth uniform Eulerian current, U, appears as a known parameter in the 
potential function, Eq. 3.1. Unfortunately, a measured pressure record gives no in- 
dication of the local current. Unlike the mean water level, the Eulerian current is a 
critical parameter that defines the propagation medium, and varying the value used 
for the current will have a considerable effect on the results of the computation. In 
fact, as pointed out by Fenton (1985a), any attempt to apply a wave theory without 
knowledge of the local current will be inaccurate, even at only first order. 

The situation is not hopeless, however. The presence of the Eulerian current 
term in the potential function draws attention to the fact that it is an important 
parameter, even if it is taken to be zero. This should encourage investigators planning 
field experiments to establish a method of estimating the local current. This may be 
as straightforward as placing a current meter nearby, or as simple as using tidal 
data to estimate the local current. Caution must be taken with the later method, 
as the accuracy of the estimate may not be consistent with the use of high order 
interpretation methods. 

It is also assumed in the LFI method that local current is depth uniform. Although 
this is unlikely to be strictly true, Kishida and Sobey (1988) demonstrated that for 


steady waves, a current profile with a realistic level of vorticity was likely to yield very 
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similar results as a depth uniform current profile for the wave generated kinematics, 
even at high order. It is expected that these results are applicable to irregular waves 
as well, and thus the assumption of depth uniform current is reasonably appropriate. 

Another complication to be considered is the direction of any local current with 
respect to the propagation direction of the waves. The U in Eq. 3.1 is the component 
of the local current in the direction of wave propagation. In the case of a single 
point measurement, no information is available to indicate the wave direction. Even 
if the current magnitude and direction are known exactly, without an estimate for the 
propagation direction of the waves, an accurate estimate for U is not available. Once 
again, further information could be used to help resolve this difficulty. Knowledge 
of the local bathymetry and the wave climate could provide this information. For 
example, waves near a beach tend to propagate almost perpendicular to the shore, or 
if there are directional measurements taken nearby available, they could be combined 
with refraction computations to provide an estimate of the local propagation direction. 
There must be some information available about both the local Eulerian current 
and the local wave propagation direction in order to interpret accurately a point 
measurement. 

The pressure records used for the examples in this chapter were generated either 
analytically or in a laboratory flume. In both cases, the Eulerian current and wave 


propagation direction are known. 


Spline of Record 


In order to avoid being restricted by the sampling rate of a particular set of 
measurements, a continuous record is computed by spline interpolation among the 
measured points. This allows the window widths and the number and location of 
samples in each window to be chosen independently of the sampling frequency of the 
data. Care must be taken, however, to be sure to include an adequate segment of 
the record in each window. A window that includes only a single measured point is 
computationally possible, but would result in misleading accuracy in the resultant 


solution. 
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Output Locations 


The spacing of the desired output locations must be chosen to determine the 
placement of the windows. The solution in each window is expected to be most 
accurate in the center, and thus a separate window is centered at each location where 
output is desired. As each window is computed independently, there is no restriction 


on the spacing of the output windows. 


Non-Dimensionalization 


The comparisons of the errors in equations of different dimensional quantities may 
result in spurious solutions that over emphasize certain equations in the formulation. 
While the familiar dimensional forms of the equations have been presented here for 
clarity, all parameters and variables are non-dimensionalized by physically identifiable 
parameters before computation. The mass density of water (p), acceleration of gravity 
(g), and the mean zero crossing frequency of the measured record (w-) are used to 


define characteristic length, time and mass scales. 


Length scale = g/w? 


Time scale = 1/w, (3.7) 
3 
Mass scale = 7 
5) 


Zz 


3.3.2 Optimization Procedure 


The primary process in the LFI method is a nonlinear optimization of a system 
typically involving up to 14 algebraic equations in 12 unknown parameters. A system 
as complex as this may have numerous local minima that result in spurious solutions. 
The best way to avoid these solutions, and to ensure efficient optimization, is to start 
with a good estimate for the unknowns in the system, and, in addition, to have a set 


of criteria for identifying spurious solutions. 


31 


Initial Estimate 


The first step in each window is to establish an initial estimate for the optimization 
procedure. Linear wave theory can be used to produce estimates for the parameters 


of correct magnitudes. The linear subsurface dynamic pressure is: 


pa = acos (kz — wt) 
cosh k(h + Zp) (3.8) 
cosh kh 


where A is the amplitude of the linear potential function. The wave number, k, and 


a= Apw 


frequency, w, are related through the linear dispersion relation: 


(w — kU)? = gk tanh kh (3.9) 


Frequency Nielsen (1986, 1989) established a method for determining the param- 
eters of a local linear approximation to waves from a pressure record. His method 
determined a local frequency and wave number that could be used to find the location 
of the water surface. A similar method is used here to determine the first estimate 
for the local frequency and wave number. Frequency of a sinusoidal signal of the form 


pa = acos(kx — wt) is available from the second derivative: 


2 
he oe (3.10) 


This approach requires an estimate for the value of the second time derivative of 
the record throughout each window. An estimate for the second derivative is directly 
available from the cubic spline of the measured points. This estimate, however, is 
very sensitive to random error in the measurements. Nielsen suggests estimating the 
value of the derivative through the use of divided differences. This solution works 
well for a smooth, sinusoidal record, but is also very sensitive to noise in the record, 
particularly in areas of small curvature; estimating a second derivative from a small 
segment of a noisy record can result in large errors. In order to accommodate the 


inevitable noise in an actual record, a different approach is taken here. 
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In each window, a third order polynomial is fit to the record in the least squares 
sense. The second derivative throughout ihe window can then be computed from this 
polynomial. By using more than four points in each window, any noise in the record 
is smoothed out by the least squares fit. This approach proved to be consistently 
reliable for artificially generated noisy records, resulting in reasonable estimates for 
the value of the second derivative throughout the window. A set of frequencies is 
computed from the estimate of the second derivative at each of the considered nodes. 
The mean of these frequencies is used as a first estimate for the local frequency of the 


record. 


Amplitude and Phase Once the frequency is known, the amplitude and spatial 
phase (kz) of a particular record can be found by rearranging the equation as a linear 


least squares problem by separating the cosine and sine components: 
Pai = acos(kzr—wt;) = bcoswt;+csinwt; (Soll) 


where a = Vb?+c? and kx = arctan(c/b). This results in the following matrix 


equation in the unknown amplitudes, 6 and c. 


coswt, sinwt, 175 
coswt, sinwt> Paz 
b : 

= (alZ)} 
Cc O 
coswt; sinwty IE ii 


? 


The system is determined if pa(t) is defined at two points. If J > 2, the system is 
over determined and can be solved in the least squared sense by common algorithms 
in numerical linear algebra libraries. This method provides a first estimate for the 


parameters: a, kz, and w. 


Refining the Linear Estimates These estimates for the parameters can be quite 


poor, as they are all dependent on the initial estimate for the second time derivative 
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of the record. In order to improve them, the estimates are refined by optimizing for 
the best least squares fit in the window. 


I 
minimize O(X) = (Pai — acos (kr — wt;))* 
xe ey (3.13) 


Xe (Gtr) 


Where P,,; is the measured dynamic pressure at time t;. The result of this optimiza- 
tion is a linear estimate for the pressure record that fits the measured record most 


closely in the window. 


Wave Number and Fourier Amplitudes Once the optimum initial frequency, 
phase and amplitude are found, the wave number is computed with the linear disper- 
sion relation (Eq. 3.9), and the first estimates for Fourier amplitudes are computed 


as follows: 


a cosh kh i 
~ pw cosh k(h + 2,) oe 
Aj =a Ay (8.115) 


The value of a is not critical and a = 0.1 was found to be satisfactory. 


Water Surface The location of the water surface at the N nodes throughout 
the window is estimated from the linear pressure response function with stretching 


(Nielsen 1989): 


Pa(tn) cosh (k (A + pa(tn)/pg)) (3.16) 


Heats pg cosh k(h + zp) 


where z, is the vertical location of the pressure gauge, and 7, is the elevation of the 


water surface node at t,. pa(t) is computed from Eq. 3.8. 
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Nonlinear Optimization 


Once there is a reasonable first estimate for all the unknown values in the potential 
function, standard nonlinear optimization routines are adequate for this system. For 
the results given, the Levenberg-Marquardt algorithm was used as implemented by 
the MATLAB Optimization Toolbox (Grace 1992). 

If the optimization routine successfully finds a minimum, the solution is checked 
to see if a clearly spurious solution is found. Spurious solutions are identified by the 


following criteria: 


e Very large or highly variable errors. 
e First order amplitude smaller than higher order amplitudes. 
e Unrealistically large or small frequency or wave number. 


e Large discontinuity between windows in the predicted water surface. 


It is unusual for the optimization routine to converge to a spurious solution. It is far 
more common for the routine not to converge at all. 

If no solution is found, or a spurious solution is found, it is necessary to revise the 
parameters of the optimization. For the next attempt, the window width is increased 
by a factor of 1.5 (1.579), and the procedure is repeated. If this is not successful, 
the window width is increased once more to twice the primary width (27). When 
increasing the window width is not successful, the order of the potential function is 
decreased until a solution is found. If none of these adjustments result in an acceptable 
solution, the window is skipped, and future analysis must be interpolated through 
that point. These adjustments are most likely to be needed in the long, flat trough 
of a shallow water wave, where the window needs to be expanded to include some 
curvature to indicate the frequency. 

Difficulties may occur in addition near zero crossings, where there is also little 
curvature in the record. Here the effects of amplitude and frequency may not be 


independent, as: 


lim asin (wt) = awt (3.17) 


t—0 
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At the limit near the zero crossing, a and w have the same effect, and the optimization 
routine can not distinguish them. Widening the window to include more of the 
surrounding record avoids this difficulty, and is generally successful in this situation 
as well as in long, flat troughs. 

Another complication can be a record that is exactly symmetric about the crest 
of a wave. In this case, the equations on either side of the crest are not independent. 
This situation is unlikely to arise in a field record, and can easily be accommodated 


by using an asymmetric distribution of points in that window. 


3.4 Theoretical Records 


To avoid complications from measurement error in the initial testing of the method, 
pressure records generated by Fourier Steady wave theory (Sobey 1989) were used. 
This approach also has the advantage of providing a solution with the complete kine- 
matics, to compare with results from the LFI method. Field or laboratory data that 
includes a full set of measured kinematics are not available. Fourier theory provides 
a near-exact solution for irrotational steady waves and can be applied at any depth 
(Rienecker and Fenton 1981; Sobey 1989). The use of analytically generated records 

‘also allows the method to be tested under a large range of conditions, by varying 
water depth, wave period, and wave height. This is useful for establishing criteria for 
choosing the solution parameters, such as window width, number of nodes in each 


window, and order of the solution. 


Shallow Water To demonstrate the optimization procedure in a window, figure 3.2 
summarizes the results from an initial estimate, before the final optimization. This 
is a window near the crest of a steep, shallow water wave generated by 18th order 
Fourier theory. The parameters of the wave are: 5m water depth, 3m wave height, 10s 
period, and zero Eulerian current, with the pressure record measured at the bottom. 
The parameters of the LFI solution are: sixth order (J = 6), and window width two 
seconds (79 = 0.27.) , centered 0.5s before the crest. 


The top plot shows the measured dynamic pressure as given by the Fourier theory 
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Figure 3.3: Final results in a window 
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generated record, and the values computed from the Bernoulli equation (Eq. 3.4), 
with the parameters of the potential function generated by the initial estimate. The 
second plot is the water surface as predicted by Fourier theory and the elevations of 
the water surface nodes generated by the initial estimate. Note that the location of 
the actual surface is given in the plot, but it is not available to help determine the 
solution. These points were all generated by the method outlined in the previous 
section, with only the pressure record as a guide. 

The third plot shows the non-dimensional errors in the objective functions: the 
Bernoulli equation and the free surface boundary conditions (Eqs. 3.2, 3.3, and 3.4). 
These are the errors that are minimized by the optimization to find the solution. If 
the solution were perfect, the error in all equations throughout the window would be 
Zero. 

The initial estimates for the dynamic pressure and the water surface have order of 
magnitude and trend agreement. The errors in the objective functions are on order 
of .03 and show a systematic pattern, particularly in the Bernoulli equation. This 
pattern, and the fact that the sharp crest of the wave has not been matched indicate 
that a better solution can be found. 

The results after the nonlinear optimization are given in Fig. 3.3. At this point the 
prediction for the dynamic pressure is essentially exact. This is virtually always the 
case, as the pressure record is available, and the solution is optimized to that record. 
The predictions for the water surface are also extremely close. This is an impressive 
achievement, as location of the water surface was found only by minimizing the errors 
in the free surface boundary conditions. The non-dimensional errors in the Bernoulli 
equation and free surface boundary conditions are on order of .001, and show no clear 
trend. The lack of trend indicates that the remaining error is random, and a good 
solution has been found. The LFI method was able to capture accurately the crest 
of a steep shallow water wave. 

Figure 3.4 shows the results of the method for the complete shallow water steady 
wave. The LFI method finds the water surface and the kinematics on the surface 
essentially exactly. While these results show the complete wave, each of the indicated 


points is in the center of a separate window, and each window was computed com- 
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Figure 3.4: LFI predictions (J = 6) and analytical kinematics at the predicted water 
surface for shallow water wave 


x Awe /g? 
© Aqw?/g? 
+ Asw?/g? 


Figure 3.5: Evolution of the solution for a shallow water wave 


pletely independently of the other windows. In this case, the parameters of the LFI 
solution are: primary window width of 2s (7 = 0.2T.), with a sixth order potential 
‘function (J = 6), seven samples on the pressure record, and seven water surface nodes 
(I = N =7), resulting in 21 equations in 16 unknowns in each window. The window 
width of two seconds is one fifth of the period of the wave, and is a reasonable length 


of time to extend the locally steady approximation. 


Evolution of the Solution Fig. 3.5 shows the evolution of the solution as the 
window is moved along the wave. The top figure shows the non-dimensional values of 
the local fundamental wave number, k, and the local fundamental frequency, w. The 
importance of the local nature of the solution is apparent in this figure, as the solution 
varies substantially from window to window. The wave number and frequency both 
increase to maximum magnitude near the crest (4 < tw, < 7), and have lower values 
in the trough region (near tw, = 3 and tw, = 9). This pattern suggests that the 


kinematics in the crest region are similar to those of a much higher frequency lower 
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order wave, and the trough kinematics similar to a lower frequency wave. 


Deep Water Figure 3.6 shows the results of the LFI method for an entire deep 
water wave. The wave was generated by 10th order Fourier wave theory, with param- 
eters: 100m water depth, 10m wave height, 10s period, and zero Eulerian current, 
with the pressure record measured 10m below the mean water level. The parameters 
of the LFI solution are: Primary window width of 1s (7 = 0.1T-), with a fourth order 
potential function (J = 4), five samples on the pressure record, and five water surface 
nodes(J = N = 5), resulting in 15 equations in 12 unknowns in each window. Once 


again, the LFI solution matches the actual solution essentially exactly. 


3.4.1 Choosing Parameters of Solution 


Unlike steady wave theory where only the order requires prior specification, this 


LFI application requires prior specification of four parameters: 
e Order of the solution (J) 
e Window width (79) 
e Number of nodes on the water surface (V) 
e Number of nodes on the pressure record (J) 


Experimentation has demonstrated that the resulting solutions are not highly sensi- 
tive to the values of these parameters; however, a reasonable solution will not result 
if these parameters are not selected judiciously. The experience acquired while devel- 
oping the LFI method on analytically generated records has provided some guidelines 
to help select appropriate values. Despite these guidelines, individual judgment must 


be used when applying the LFI method to a particular measured record. 


Number of Nodes on the Water Surface and Pressure Record Mathemat- 
ically, the system of equations is specified if there are at least as many equations as 


unknown parameters. In this case, /+ N > 4+ J. While meeting this criterion is 
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Figure 3.6: LFI predictions (J = 4) and analytical kinematics at the predicted water 
surface for deep water wave 
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the minimum mathematical requirement, there are other factors that must be taken 
into account to assure a reasonable physical solution. The local water surface is rep- 
resented by N nodes in each window. Defining the surface with N = J+ 1 points 
achieves a representation for the water surface of the same order as the potential 
function being used. Care must be taken with very low order solutions; for example, 
2 points would define the water surface to first order, but would provide no indication 
of the curvature of the surface. As a rule of thumb, at least three points should be 
used, regardless of order. 

Defining the water surface at J +1 points provides 2(J + 1) equations (f? and 
f*, at the J+1 points). To keep the order of approximation consistent, the Bernoulli 
equation is also applied at J+1 points within each window. Using J+1 water surface 
and pressure record nodes overspecifies the solution at all orders. This arrangement 
of equations proved effective for all the examples given in this chapter on analytically 
derived records. 

While that approach was effective on these few examples, it is important that the 
final solution matches the measured record closely. It’s possible that this suggested 
arrangement of points would allow the optimization routine to be biased towards 
the more numerous free surface boundary conditions, giving a solution that does not 
match the pressure record well. The Bernoulli equation could be applied at more 
points on the pressure record than the number of nodes defining the water surface. 
Increasing the number of points at which the Bernoulli equation is applied will shift 
the emphasis of the optimization to the measured record. Additional points on the 
measured record can also be useful for accommodating measurement noise that may 


be present in field or laboratory records. 


Order of Solution Similarly to high order steady wave theory, the order chosen 
for the solution is influenced by a number of factors including the height of the waves, 
the depth of the water, and the accuracy desired. As with steady wave theory, higher 
order results in greater accuracy at the expense of computational simplicity, and is 
necessary for larger waves and for shallow water. The following examples will help to 


provide guidelines for the order chosen. 
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Figure 3.7: Crest of a shallow water wave at order 1 
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Figure 3.8: Crest of a shallow water wave at order 2 
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Figure 3.9: Crest of a shallow water wave at order 3 
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Figure 3.10: Crest of a shallow water wave at order 4 
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Figure 3.12: Crest of a shallow water wave at order 6 
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Figures 3.7 through 3.12 show the crest of the same shallow water wave shown in 
Fig. 3.4, computed to orders J = 1 through 6. At first order, the non-dimensional 
errors in the Boundary conditions and the Bernoulli equation are of order 10~°. These 
are small errors, but the sharp curvature at the crest has not been captured. As 
the order is increased, the magnitude of the errors increases. The increase in the 
magnitude of the errors is due to the increase in the number of equations included. 
There are three additional equations included in the optimization for each increase 
in order (the two FSBCs at each additional water surface node, and one additional 
pressure record node), but there is only one more parameter in the solution vector. 
In order to best satisfy this system, the optimization finds a solution that results in 
slightly more error in each equation. The magnitude of the errors does increase, but 
the solution slowly converges to very precisely match the sharp crest at sixth order. 
An asymmetric distribution of points was used in this window to accommodate the 
symmetry about the crest. 

Figures 3.13 through 3.16 show the crest of the same deep water wave as Fig. 3.6, 
computed to orders 1 through 4. Even at first order, Fig. 3.13, the solution is very 
good. Deep water waves generally do not require a very high order solution, linear 
wave theory often being reasonably adequate in these conditions. It’s important to 
keep in mind that, although this solution is first order, it is still nonlinear, having 
found a minimum in the errors of the full nonlinear governing equations, and the 
frequency and wave number are free to vary, not being bound by the linear dispersion 
relationship. The local nature of the solution would allow it to change with time, 
accommodating an irregular profile at low order better than an equally low order 
global solution. 

In this case, first order provides an acceptable solution; however, the water surface 
is more accurately matched as the order is increased, and the solution converged well 
at higher order, so there is little penalty in using up to fourth order. For an irregular 
record, higher order is more likely to be successful in matching the irregularity in the 
record. As the examples show, deep water waves can be well represented at low order, 


while higher order in necessary in shallow water. 
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Figure 3.14: Crest of a deep water wave at order 2 
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Figure 3.16: Crest of a deep water wave at order 4 
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Window Width Choice of window width is a balance between including sufficient 
curvature in the record to identify the local wave response while minimizing the extent 
of the locally steady approximation. Another factor to be considered is the sampling 
rate of the data. The LFI method employs a cubic spline for interpolation among 
the measured points. This approach allows the selection window widths and sampling 
locations without being restricted by the data sampling locations. While this freedom 
is very useful, it can be deceiving, as the data are not truly continuous. As the data 
are not continuous, it is important to be sure that the window width chosen includes 
sufficient measured data points to justify the order being used. Unfortunately, field 
data are often sampled at a frequency of as low as 1Hz. In this case, a window width 
of one tenth the mean zero crossing might be on order of 1 second. This would provide 
a maximum of two actual data points in each window. It may not be appropriate 
to attempt a high order solution with such little hard data. Ideally, the LFI method 
should be used with data sampled at a higher rate. If that isn’t possible, wider 
windows, and perhaps lower order solutions, should be used. 

In the case of the theoretical records used in the previous section, the data sam- 
pling rate could be arbitrarily high. Without being restricted by measurement fre- 
quency, some guidelines for window width have been determined. A primary window 
' width of one tenth the zero crossing period provides a good starting point. Using 
a smaller window often did not allow convergence of the solution. It was necessary 
to occasionally increase window width at some locations on the record, as discussed 
in section 3.3.2. When a primary window width is determined that works for most 
of the record, while needing to be increased at a only a few locations, the optimum 
width has been found. 

In the case of the shallow water wave discussed in the previous section, a window 
width of 0.17. was successful for up to a fourth order solution. Unfortunately, the 
fourth order solution did not fully capture the sharp crest of the wave. When a higher 
order solution was attempted, the optimization would not converge with the given 
window width. A wider window had to be used in order to obtain a solution of high 
enough order to capture the sharp crest of the wave, and thus a window width of 0.27, 


was used in that example. In general, it has been necessary to use wider windows for 
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higher order solutions. 


3.5 Laboratory Measurements 


The following results use laboratory data collected by Murray Townsend and 
John Fenton at the Australian Maritime Engineering Cooperative Research Center at 
Monash University, Melbourne, Australia in June, 1996. The dynamic pressure, fluid 
velocities and water surface were measured at the same horizontal location along the 
flume; the pressure at a variety of elevations, and the horizontal and vertical velocities 
at an elevation of -0.8 meters, with the origin at the mean water level. The still water 
level was 1.55 meters, with the data sampled at 60 Hz. 

The wave flume is 52m long by 2.2m wide with two working sections 4m and 2.5m 
long connected by a ramp. A false floor had been built into the shallower working 
section to bring the depth to 1.55m. At one end of the flume is a hydraulically oper- 
ated bottom hinged wave paddle and at the other a beach that absorbs a minimum 
of 96% of wave energy across the range of operating frequencies. The measurements 
were taken approximately 22m from the beach end of the flume. The beach length is 
6m. 

The wave surface was measured with a capacitance wave probe with typical cali- 
brations providing a correlation coefficient R? greater than 0.999. The pressure was 
measured with Druck PDCR35/D transducers located in the flume near the center 
of a long (2.5m long by 2m high) plywood board aligned with the direction of wave 
propagation. The measurement face of the probe was flush with the board surface 
to eliminate local flow separation. The R? values from calibrations were in all cases 
above 0.999. Fluid velocities were measured with a factory calibrated Alec electronics 
ACM250-A electromagnetic current meter which is capable of velocity measurements 
in two dimensions. 

One of the difficulties encountered with pressure measurements is the question 
of what has actually been measured. If the transducer itself is absolutely accurate, 
the pressure recorded will be the actual pressure at the location of the transducer. 


Unfortunately, the presence of the instrument is likely to alter the flow in its vicinity, 
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creating substantial differences between the pressure at the transducer, and the pres- 
sure that would exist if the transducer were not there to disturb the flow. Another 
source of error is the frequency response of the transducer (Raichlen et al. 1990; 
Ippen and Raichlen 1957). The frequency response of pressure transducers is not flat, 
and as a result, the recorded signal could be quite different from the actual pressure. 

In the experimental setup described above, the pressure transducers were mounted 
on a large plywood panel oriented in line with the flume. The measurement face of 
the instrument was flush with the surface of the plywood to minimize dynamic effects 
near the gauge. The panel was large enough for the boundary layer to be fully 
developed near the surface of the plywood, to prevent edge effects from the edge 
of the plywood affecting the measurements. This arrangement is expected to have 
resulted in minimal dynamic effects on the measured pressure. 

There is no information available about the frequency response of the pressure 
transducers used for this experiment. However, a spectral analysis of the records can 
help identify any potential problems. The top plot of figure 3.17 gives the measured 
water surface and subsurface pressure for a short segment of a record. There are 
some clear higher frequency fluctuations in the water surface that do not appear in 
the pressure record. In this particular segment, there is a sharp secondary crest in 
the first trough (near the 2s point) in the water surface record. 

This loss of high frequency information in the pressure records is confirmed by 
an examination of the discrete Fourier transform of the records. The second plot in 
figure 3.17 is the smoothed variance spectra of the two records from which the above 
segments were taken. The water surface record has a great deal more variance at 
the higher frequencies. Note that the spectrum of the pressure record has decayed to 
almost zero at w = 4s~! point, while there is still considerable variance in the water 
surface at that frequency. 

Linear wave theory predicts that the motion of high frequency waves decays with 
depth much faster than low frequencies. In this example, the non-dimensional water 
depth at the peak frequency is w*h/g ~ 1. This is generally considered intermediate 
depth water, and the decay with depth is expected to be moderate. At the frequency 


where there is essentially no energy in the pressure record, the non-dimensional water 
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Figure 3.17: Water surface, pressure record, and variance spectra of records from 
flume experiment 
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depth is w?h/g =~ 2.5, which is considered to be deep water. In deep water, the 
action of the waves decays very quickly with depth, so that little energy is expected 
to remain at half the water depth, the depth at which the pressure measurements 
were taken. 

The third plot in figure 3.17 is the spectrum of the subsurface pressure record. 
The solid line was computed from the measured record, and the dashed line was 
computed from the water surface record, using the linear pressure response factor: 


cosh k(h + z,) 


ap(w) = a,(w) cosh( Eh) (3.18) 


where a,(w) and a,(w) are the magnitude of the pressure and water surface amplitude 
spectra at a given frequency, and the wave number, k, is computed from the linear 


dispersion relationship: 
(w — kU)? = gk tanh(kh) (3.19) 


The predicted spectrum matches the measured spectrum very closely. This indicates 
that the loss of high frequency information in the subsurface pressure record is the 
result of the predicted decay with depth of the energy at the higher frequencies, rather 
than a result of the frequency response of the pressure gauge used. 

The final plot in figure 3.17 is the spectrum of the water surface record. The 
solid line was computed from the measured water surface, and the dashed line was 
computed from the pressure record, also using the linear pressure response factor. The 
predicted spectrum of the water surface matches the measured spectrum well near 
the peak frequency, but strongly over-predicts the high frequencies. If this method 
were used to predict the water surface, a somewhat arbitrary high frequency cut off 
would have to be established. 

An examination of the top plot in figure 3.17 reveals what appears to be two 
dominant free modes in the water surface. There is a large low frequency mode 
with a period of approximately 2s, and a higher frequency mode with a period of 
approximately 1s. If the LFI method were applied to the water surface record (Sobey 
1992), the moving widow could capture a different dominant mode in each window. 


Unfortunately, the higher frequency mode has decayed with depth too much to be 
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Figure 3.18: Measured and predicted water surface near the crest of a sharp wave in 
the flume 


clear in the pressure record; applying the LFI method to the pressure record results 
in predictions that capture the dominant low frequency mode, and miss the higher 
frequency mode. If the higher frequency information is not in the local segment of 
the record, it is not possible to reconstruct it. 

In the previous section, the LFI method was applied to pressure records generated 
by Fourier steady wave theory. These records had only a single free mode, and thus 
the LFI method was able to identify that mode. In addition, in shallow water, there 
is little decay with depth of the wave action, and the LFI method could be expected 
to be effective, as it was on the Fourier record. In deep water, the method can be 
effective if the pressure gauges are located near enough to the surface. For the deep 
water steady wave presented previously (figure 3.6) the pressure record was taken 
from a depth of one wave height below the water surface, one tenth of the total 
depth. At this depth, there is adequate information to identify the action at a wide 
range of frequencies. 

Given the limitations of the data, a segment of the record was chosen in which 


there did not appear to be a substantial high frequency component to the water 
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Figure 3.19: Segment of record from flume measurements 
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surface record, and thus the effect of the high frequency decay should cause few 
problems, and the LFI method can be applied. 

While this particular segment was chosen to minimize the effects of the high 
frequency decay, it was not possible to eliminate it entirely. Figure 3.18 shows the 
measured water surface and the LFI predictions of the water surface in three windows 
in the vicinity of a sharp crest in the record. The solution parameters are: J = 3, 
™ = 0.3T., and J = N = 4. Both the pressure and the velocities were measured 
at an elevation of 0.8m below the MWL, with the mean water depth 1.55m. The 
LFI method located the water surface fairly accurately, although it clearly missed 
the very sharp peak of the crest. Varying the solution parameters did not improve 
the accuracy at the crest. In fact, a smaller window width resulted in a discontinuity 
in the water surface between windows. The inability of the LFI method to capture 
the sharp crest is likely due to the missing high frequency information. If the high 
frequency part of the signal is missing from the examined pressure record, there is no 
way to recapture that information. In order for the method to be more effective in 
this situation, the pressure would have to be measured at a depth closer to the water 
surface. 

Figure 3.19 shows the results of the LFI method applied to the entire segment of 
the record from which the crest was taken. The water surface is captured fairly well, 


as are both the vertical and horizontal velocities. 


3.6 Discussion 


The given results demonstrate the potential of the LFI method in the interpre- 
tation of submerged pressure traces in a variety of conditions. In the case of steady 
waves, the LFI method accurately computed the detail of the wave, using only data 
from a small window in time. In particular, the method was able to capture the 
pronounced sharp crest of a steep, shallow water wave. 

The method did not perform as well on laboratory records, failing to capture 
some of high frequency detail in the water surface. This is due to the limitations of 


subsurface pressure records, where much of the high frequency information is missing 


of 


from the record. The method is likely to perform better in shallow water, or with 
records that are measured closer to the water surface. Despite this limitation, the 
method was able to capture much of the detail of a irregular laboratory record, and 
provide fairly accurate estimations of the local kinematics. 

The analysis of regular waves provides guidelines for the parameters to be used 
in the analysis of irregular waves. In shallow water, higher order solutions and wider 
windows must be used than in deep water. Window widths of one fifth of the zero 
crossing period and a sixth order potential function are adequate for the shallowest 
waves, and window widths as small as one tenth of the zero crossing period and a 
third order potential function are adequate for deep water. 

The laboratory results indicate that the method requires good precision and care in 
the measurements. Any high order method demands very accurate data, but this need 
is exacerbated by the ill posed problem of determining the near surface kinematics 
from a subsurface record. While not particularly sensitive to random noise in the 
record, the decay with depth of the high frequency information makes it very difficult 
to capture the high frequency modes. Fundamentally, the LFI method is designed to 
capture the dominant free mode in each given window. As the higher frequency modes 
decay faster with depth, and if the measurements are taken far below the surface, the 
dominant mode will always be one of lower frequency modes. This difficulty would 
be exacerbated by any limitations in the frequency response of the gauges. 

The computer resources required for the method are substantial, but not pro- 
hibitive. As computers continue to get faster, computation time is not the considera- 
tion that it once was. The method was developed and all computation done using the 
MATLAB computational package. MATLAB is an interpreted language that provides 
a very easy to use interface to a robust and complete library of computational and 
visualization routines, allowing for rapid development of new methods. Being an in- 
terpreted language, the resulting code is not as fast as it might be if the routines were 
programmed in a compiled language, such as Fortran or C. The computational speed 
is also affected by the degree to which the program pauses to provide visual output. 
Perhaps a better measure of the computational intensity of the method is a count of 


the total number floating point operations (flops) used to compute the solution. 
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The examples in this dissertation were all computed on a Intel 90MHz Pentium 
processor based PC with 32MB of RAM, running the Linux operating system. As an 
order of magnitude estimate, the shortest computation time for an example in this 
chapter is 4.8 minutes and 7.8 x 10° flops for the complete laboratory record presented 
in figure 3.19. The longest computation time required was 71 minutes and 582 x 10° 
flops for the shallow water wave presented in figure 3.4. There are a number of 
reasons for this large variation in computation time. The first is simply the number 
of window solutions computed. As each window solution is computed separately, 
the computation time increases linearly with the number of windows computed. If 
computational time is a concern, this can be taken into account when choosing the 
output spacing. The other reason that the shallow water wave takes much longer 
to compute is that there are a number of windows that did not converge with the 
first set of computational parameters. The optimization is run for many iterations 
to ensure that it won’t converge. As the parameters are varied, the computation is 
repeated. This process takes a great deal of time. 

With further development, it may be possible to determine a set of criteria for the 
computational parameters by examining the segment of the record in a given window. 
This would be far more efficient than simply attempting a solution with a variety of 


values until convergence is achieved. 
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Chapter 4 


Three Dimensional LFI Theory 


The previous chapters were concerned with the determination of wave kinematics 
from the analysis of a time series of a single physical quantity at a single location. By 
assuming that the flow is two dimensional, a reasonable approximation of the wave 
kinematics can be found. Unfortunately, this analysis does not give any information 
about the directionality of the sea state. There are some processes in which the wave 
directions are directly important, such as sediment transport. Even in situations 
where the wave directions are not directly important, it has been shown that omitting 
_ the directional nature of the sea results in substantial errors in the prediction of scalar 
quantities, such as the maximum velocities and accelerations in a measured wave crest, 
or the resultant forcing on structures (Dean 1977; Forristall et al. 1978). 

In order to capture the directional nature of the sea, an array of instruments must 
be used. The result is a set of time series of a single physical quantity at a number 
of different locations, or a set of different physical quantities at the same or different 
locations. This chapter outlines a method for determining the directional kinematics 
of irregular seas that can be adapted to accommodate virtually any combination of 


such time series. 
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Figure 4.1: Coordinate system used for 3-d method 


4.1 Three Dimensional Seas 


The development of the two dimensional LF] method was closely anchored to 
the very complete understanding of two dimensional steady waves. In contrast, the 
understanding of three dimensional wave fields is not nearly as complete. Much of 
the literature on three dimensional seas attempts to describe the motion through the 
use of a directional energy spectrum. Far less attention has been directed to the 


determination of the detailed kinematics of directional seas. 


4.2 Problem Formulation 


The governing equations for three dimensional gravity waves are a straightforward 
extension of those in two dimensions to include the third dimension. The flow is taken 
to be irrotational and incompressible, and thus the kinematics can be represented by 
a potential function, d(x, y, z,t), in a Cartesian coordinate system (Fig. 4.1), where: 

_ Of _ Od (aka) 


u and v are the horizontal velocities in the z and y directions, and w is the vertical 


velocity. 
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Field Equation The field equation is mass conservation for irrotational flow, rep- 


resented by the Laplace equation: 


Od Oo O76 
Vv? = a 
Me Ox? i Oy? x Oz? 


0 (4.2) 


Boundary Conditions The boundary conditions are the bottom boundary condi- 
tion (BBC) for a horizontal bed, 


as 
=> = 


W 


0 at z=—h (4.3) 
the dynamic free surface boundary condition (DFSBC), 


1 
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where 77 is the elevation of the water surface. 

As with the two dimensional formulation, the kinematic free surface boundary 
condition requires the gradients of the water surface (22, on and a): Knowledge 
- of these gradients is likely to be limited. In the case of an array of water surface 
measurements, estimates for the gradients could be computed by interpolating among 
the measured elevations, but this would result in a low order estimate, and compound 
the inevitable error in the measurements. In the case of subsurface measurements, 
there are no data as to the location or the gradients of the water surface. 


To eliminate these gradients, a modified kinematic free surface boundary condition 


is defined by differentiating the DFSBC following the motion (Longuet-Higgins 1962), 


as was done in two dimensions in Section 2.1. 


MKFSBC = —g(KFSBC) + =. (DFSBC) = 
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This form does not include the gradients of the water surface, and all the terms are 
defined by the potential function. Applying both the modified kinematic free surface 
boundary condition and the dynamic free surface boundary condition completes the 


formulation without the need for knowledge of the gradients of the water surface. 


Observational Equations The field equation and the boundary conditions de- 
scribe a free surface potential flow. In steady two dimensional wave theory, a wave 
height and periodicity in space and time are specified to complete the formulation. 
In the two dimensional LFI theory (Chapter 2), a form for the potential function is 
chosen, with the parameters determined to fit the free surface boundary conditions 
and the measured record in a small window in time. In three dimensions, the pro- 
cess is essentially the same. A three dimensional form for the potential function will 
be presented, with parameters that are found to fit the measured records and the 
boundary conditions. 

In order to define a solution that fits the measured record, observational equations 
are established. These equations are defined to make use of the particular quantities 
that have been measured. In the case of an array of water surface measurements, they 
are the free surface boundary conditions, applied at the measured elevations and hor- 
izontal locations at a number of points in time throughout the window (Chapter 5). 
In the case of an array of pressure measurements, they are the Bernoulli equation, 


applied at the elevation and horizontal locations of the measurements, and also at 
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a number of points in time within the window (Chapter 6). The method could be 
adapted to virtually any combination of measured physical quantities by establishing 
appropriate observational equations. In addition to the aforementioned water surface 
and pressure measurements, these could include water surface gradient or accelera- 
tions measured by wave buoys, subsurface velocity measurements, or any combination 


of these. 


4.3 Formulation of the Solution in each Window 


4.3.1 Background 


Two dimensional steady wave theory provided the inspiration for the development 
of the two dimensional LFI method (Chapter 2). Unfortunately, the literature does 
not provide as solid a basis for nonlinear interpretations of directional seas. There 
has, however, been some work that can be used as a basis for a directional LFI 
method. In an early attempt to explore the nonlinear nature of directional seas, 
Longuet-Higgins (1962) computed the interaction of two intersecting steady waves 
in deep water through the use of a double perturbation expansion in the steepness 
_of the waves up to third order. The result was a potential function that contained 
terms representing the higher order interaction between the phases of the intersecting 


waves: 
bo” = fi(z, 51) + fo(z, So) 
po) = fs(z, 51 + 52) + faz, 51 — So) 
$°) = f5(z, 251 + S2) + fe(z, 2:51 — S2) + fr(z, Si + 252) + f(z, 51 — 252) 
Sn = (Kn -X — Wnt + An) 


(4.7) 


where n = [1,2], d'”) is the mth order potential function, S; and S> are the phase 
functions of the two waves, k,, is the vector wave number of wave n, x is the horizontal 
position vector, w, and a, are the angular frequency and initial phase of the nth 
wave. The functions, f;, were determined algebraically by expanding the free surface 


boundary conditions in a Taylor series about the mean water level, and solving for 
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the potential function at each order. Nonlinear frequency modulation was not taken 
into account, so the frequencies and wave numbers of each wave independently satisfy 
the linear dispersion relation for water of infinite depth: 
w, = glk,| (4.8) 

fi and fz are independent of each other, and are the familiar linear wave solution: 

ja S erie sin $1 i mre sin So (4.9) 
The higher order terms are all functions of both waves, and thus include the interac- 
tion of the two waves. 

A number of investigators subsequently expanded upon this work. Hsu and Chen 
(1992) presented a detailed examination of Longuet-Higgins (1962), pointing out diffi- 
culties that arise from the assumption of the linear dispersion relation. They presented 
a more mature analysis, including higher order modulation of the wave frequencies, 
and higher order self interactions of the individual waves. This resulted in a com- 
plete theory up to third order for two intersecting waves in deep water. Hsu and 
Chen also proposed a systematic ordering in the phase relationships generated by the 
interactions of two steady waves to arbitrary order. 

Expanding upon this work, Ohyama, Jeng and Hsu (1995a) extended the per- 
turbation expansion method in a number of ways. The most recent version of the 
method accommodates any number of waves, allows for water of finite depth, and is 
accurate to fourth order. This last method can compute the water surface and full 
kinematics of a highly irregular sea, as produced by a large number of intersecting 
waves. A more detailed discussion of the method is given in Appendix A. 

Ohyama, Jeng and Hsu’s work suggests a form for the potential function that can 
be used for a local method for the recreation of three dimensional kinematics from 
arrays of wave measurement devices, including water surface arrays, pressure arrays, 
directional current meters, or any combination of these. 


The general form for the potential function representing N intersecting steady 
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waves at order J is: 


P(z,y,2,t) = Urx st Uyy ate 


J=i1| J=TNTt lim| 
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where U, and U, are the components of the depth uniform Eulerian current, k,,, and 
ky, are the vector wave number components of the nth wave in the z and y directions. 

There is a separate summation for each wave considered. 
As the notation in Eq. 4.10 is quite confusing, examples with two and three waves 


are given below. With two waves (N = 2): 


P(2,Y,2,t) = Ure + Uyy + 
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And with three waves (N = 3): 
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Sn = (Kent + ky ny — Wnt + On) fe = [1,23] 


(t1k2.i: + takzin + t3kz,ig)? + (21k yi: + 22h yin + takty,is )? 


This form for the potential function exactly satisfies mass conservation, in the form 
of the Laplace equation, and the bottom boundary condition for a locally horizontal 
bottom. It allows for high order representation of each of the steady waves, as well 
as the interaction of each wave with every other wave. The balance of the solution is 
specified by the full three dimensional free surface boundary conditions, Eqs. 4.4 and 


4.6. 


4.3.2 Dynamics 


The dynamics are available through the Bernoulli equation in three dimensions: 


ag 
ot 


1 =e 
shag ts aa a =1() (4.13) 


The dynamic pressure is the difference between the total and hydrostatic pressure 


(Pa = P — Ph)- 


The Bernoulli Constant 


In Chapter 2, an explicit and exact expression for the Bernoulli constant, B, is 
given for steady two dimensional waves. (Eqs. 2.14 and 2.15). A similar expression 
can be established for unsteady three dimensional waves. 

The Bernoulli equation is applied at the bed: 


Od, 1 i >) 
a 7 CO) eos 2 (4.14) 
where the subscript, 6 indicates the value at the bed. 

The two dimensional approach is based on analysis first presented by Longuet- 
Higgins (1975). In that work, the analysis was applied to steady waves. With steady 
waves, the average over either time or space of the pressure at the bed is the hydro- 


static pressure. With unsteady waves, the pressure at the bed averaged over space 


at any given time, or averaged over time at any given location, is not necessarily the 
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hydrostatic pressure. If, however, an average is taken over both time and horizontal 


space, the result must be the hydrostatic pressure. 


b = pg(7] — 2) (4.15) 


3 


with the double overline indicating an average over time and horizontal space. When 
the Bernoulli equation at the bed is averaged over time and horizontal space, the 
average time gradient of the potential function is zero, resulting in a simple expression 


for the Bernoulli constant: 
(u2+ v7) +97 (4.16) 


For the potential function given by Eq. 4.10, and z = 0 at the mean water level 


(7), B becomes 
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giving a complete expression for B as a function of the parameters of the potential 


function. 


4.4 A Local Two-Intersecting-Wave Theory 


While a large number of intersecting waves could capture a sea state of virtually 
any complexity, it would be difficult to distinguish between the effects of each in- 
dividual wave. It is important to remember that the familiar directional spectrum 
description of a sea state is an integral description. While there are many different 
frequency and direction modes represented in the spectrum, they do not necessarily 
all play a significant rule in the kinematics the entire time. In fact, when observing 
irregular seas, it is often the case that at any given time, there appears to be a single 
dominant wave of a particular frequency and direction. Over time, there is a series of 
such dominant waves, each of a different frequency and direction. The integral effect 


of this process is a broad directional spectrum, but if time is separated into individual 
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small windows, in each window there may be only one or two free waves dominating 
the motion. 

Taking advantage of the local nature of the approach, only one or two intersect- 
ing free waves are considered in each window. While unlikely to capture the full 
complexity of a broadly directional sea, this approach should be able to capture the 
dominant modes in each window. The direction, frequency, and amplitude of the 
dominant modes will vary substantially from window to window as they do in the 
two dimensional approach, having an integral effect that includes a large variation in 
directionality. 

The two wave method is expected to be particularly effective for the case of 
standing waves or short crested waves, as would be found near a reflecting surface 
when the incident wave field is almost unidirectional, as it often is in shallow water. 


Expanding Eq. 4.11 to fourth order, the potential function for two intersecting 


waves is: 
20 
P(x, Y, B..h) =U,2+ Uyy + Sy A; (KA, Kz) sino; (4.18) 
cl 
where: 
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At third order: 


o7 = (351) og = (352) 
Og = (25; + 52) G19 = (25; — So) 
011 = (Sy + 252) 012 = (5) = 2S2) 


13 = (45)) O14 = (452) 

O15 = (35; + $9) Cig = (35; — $2) 
O17 = (25 + 252) o1g = (25, — 22) 
O19 = (5; + 352) 729 = (51 — 3S) 


01,02,03,04,07,08, 013, and O14 are self-wave interactions. The balance of the o; are 
Wwave-wave interactions. 

Because the phases are arguments of the sine in the potential function, a sim- 
ple expansion of Eq. 4.11 results in redundant terms that have been removed. For 


example, if there are two components: 
By sin(S;—S2) and Bposin(S2—5;) (4.19) 
Because sine is an odd function: 
By sin(S2 — $,) = —Bz sin($; — S2) (4.20) 
the two components can be combined into a single component: 
(B, — B2)sin(.S; — Sz) = A;sin( Si — S2) (4.21) 


If both components were included, the effects of B; and B2 would be indistinguishable 
in the optimization; they must be combined into the single coefficient, A;. 
This form for the potential function results in 28 unknowns at fourth order (20A; 


plus kz, ky, w, and a for each of the two intersecting waves) that must be found 
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for each window in time. At least 28 independent equations must be defined by 
applying a combination of each of the free surface boundary conditions (Eq. 4.4 and 
4.6) and the observational equations at a number of nodes throughout each window. 
The specific combination of these equations will be determined by the layout of the 


instruments. 


4.4.1 Kinematics 


In order to apply the free surface boundary conditions, the full kinematics must 
be known. While the kinematics are completely specified by the potential function, 


Eq. 4.18, some of the algebra may not be obvious. The full expressions are provided 


here. 
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where: 


sees 
CORBI ae 


cosh h;h 
oa ey , Sul AG((lo <r 2) 
C(KA, K;z) = ~ aoa 
Paes 00; 2 4 Oo; g 
a Ox Oy 
Oo; 
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00; Oo; 
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The solution details for a given set of measurements are problem specific. These 
will be provided for an array of water surface measurements in Chapter 5 and for an 


array of pressure measurements in Chapter 6. 
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Chapter 5 
Array of Water Surface Traces 


Water surface measurements are difficult to obtain and not very commonly used 
for field measurements, but they are routinely used in the laboratory. Surface piercing 
wave gauges are the most commonly and easily deployed of the methods for measuring 
waves in the lab. Unfortunately, even large arrays of such wave gauges still only give 
information about the fluctuations of the water surface. The underlying kinematics 
must still be predicted with a wave theory. The following is a method for determining 
the kinematics of waves in the region of an array of water surface measurements. It 
is an extension of the LFI method presented in the previous chapters, expanded! to 
determine the directional structure of the wave field. The method is fully nonlinear, 
and results in a complete prediction for the full kinematics in the vicinity of an array 


of water surface measurements, throughout the depth of the water column. 


5.1 Formulation of Solution 


As described in more detail in Chapter 4, the flow is assumed to be irrotational and 
incompressible, with a potential function that represents either a single directional 
wave, or two separate intersecting waves. 

When the measurement is taken in a location that is far from reflecting surfaces, 
it can be effective and straightforward to assume that the wave field can locally 


be defined as a segment of a single steady wave. This is quite similar to the LFI 
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method for a point water surface record (Sobey 1992), but with the direction of the 
wave determined. The varying directional trend of the sea state is accommodated 
by determining the direction in each window separately, allowing for a fluctuation of 
the wave direction with time, from window to window. In this case, the potential 
function is Eq. 4.10 reduced to a single wave: 
cosh 7A (h + z) 
cosh jh 


IX = ene 


where U, and U, are the components of the known depth uniform Eulerian current 


J 
(z,z,t)=U,2 + Uyy + Se A; 


j=l 


sinj(kpz +kyy—wt+a) (5.1) 


in the z and y directions, h is the mean water depth, J is the truncation order of the 
Fourier series, A; are the Fourier coefficients, w is the local fundamental frequency, 
k, and ky, are the components of the local wave number in the z and y directions, 
and Kk is the magnitude of the local wave number. 

When an array of water surface gauges is placed near a reflecting surface, such 
as a sea wall, the resulting sea state is likely to contain simultaneous components in 
different directions, such as in a standing wave or short crested sea. This effect can 
not be captured with a potential function representing a single progressive wave. It is 
~ possible, however, to capture this type of sea with a potential function representing 
two intersecting waves. The potential function for two intersecting waves is Eq. 4.18, 


which is repeated here: 
20 
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a1 = (51) 
At second order: 
oe, = (2Sn) 
65 = (S; + $2) 
At third order: 
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And at fourth order: 
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It might be possible to include a larger number of intersecting free waves in order 


to capture a broadly directional sea, but it would introduce additional complication 


in distinguishing the effects of the individual waves, and will not be considered here. 


Both of these potential functions satisfy the field equation (Laplace, Eq. 4.2) and 


the bottom boundary condition (Eq. 4.3) exactly. The balance of the solution is 


determined by the free surface boundary conditions: the modified kinematic free 
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surface boundary condition (f*), 
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with the Bernoulli constant (B) defined as: 
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The problem of determining the kinematics of irregular waves from a set of mea- 
sured water surface traces is a mathematically better posed problem than interpreting 
a subsurface record. The flow is governed by the elliptic Laplace equation (Eq. 4.2), 
so that the solution is determined by the boundaries. While the complete boundaries 
of the solution domain are not known, the boundary conditions and location of the 
boundary are known at both the top and bottom of the solution domain. The bottom 
boundary condition is well defined, and the location of the water surface is measured 
at a few locations in space and many points in time, allowing the direct application of 
the free surface boundary conditions. This is in contrast to working with subsurface 
records, in which the location of the free surface must be determined in order to apply 
the free surface boundary conditions. The need for horizontal boundary conditions is 
eliminated by the assumed periodicity of the chosen potential function. However the 
fundamental wavelength(s) and period(s) must be found as part of the solution. 

When working with a point measurement, the fundamental frequency is fairly well 
defined by the time evolution in the window chosen, as long as the window is wide 
enough. There is no direct information available about the spatial evolution of the 


signal, however, so the wave number is determined only through the application of 
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the free surface boundary conditions. An array of measurements provides information 
about the spatial evolution of the signal, helping to better define the fundamental 
wave number. This usually results in faster and more robust convergence of the 


optimization. 


5.2 Formulation of the Optimization 


The formulation of the optimization for the LFI method as applied to the in- 
terpretation of an array of water surface measurements has a great deal in common 
with the formulation for a subsurface pressure record (chapter 3) and a point water 
surface measurement (Sobey 1992). Less detail is presented here than in the previous 
chapter, but the framework will be presented, with an emphasis on the additional 


information necessary for applying the method to an array of measurements. 


Observational Equations The governing equations presented in the previous sec- 
tion represent a free surface potential flow, with one or two components propagating 
in an arbitrary direction. The observational equations are the equations in the system 
that force the solution to fit the given measured record. As the location of the water 
surface has been measured, these are the free surface boundary conditions (Eqs. 5.3 
and 5.4), applied at the horizontal location of each of the nodes in the array. Sufficient 
independent equations are defined by applying the boundary conditions at a number 
of times along the measured records, within the window in time considered (Fig. 5.1). 
The solution is the set of parameters in the potential function that result in the least 
error in the FSBCs. 

In order to specify the solution, there must be at least as many independent 
equations as there are unknown parameters in the system of equations. The free 


K 


m,n 


surface boundary conditions ( and f?_) are applied at each of the N measured 
locations and at M time samples in the window, resulting in 2MN independent 
equations. In the single wave case, there are 4 + J unknown parameters sought in 
Kq. 5.1 (kz, ky, w, a, and A,...A,z) in each window, so that if 2MN > 4+ J, the 


solution is specified. In the two wave case, there are 2)> 7 +8 unknowns in Eq. 5.2 


Ue 


t 


Figure 5.1: Schematic of system of equations in a window 


(2k, 2ky, 2w, 2a, and A,...Azy: 10 at first order, 14 at second order, 20 at third 
order, and 28 at fourth order). The solution is specified when 2MN > 2507 +8. 


This system results in the following least squares optimization. 


D 2 
minimize O(X ass SE ( 1G a Ma Uncen SPS Oa Ontbnenta) (GS) 


m=1 n=1 


where: 
X= (ke Sp Ds Qa, vAre D009 Aj) 
for the one wave case, and: 
X = (koi, ky1, 1, 01, ke,2, ky,2,W2, 02, Ar,--. , Az) 


for the two wave case. z, and y, are the coordinates of the nth gauge, and mn is 
the measured elevation of the water surface at time t,, and gauge n. 

As with the analysis of a pressure record, overspecification can be helpful in accom- 
modating the measurement errors in field records. More than the minimum number 
of time samples in the window may be required to define the shape of the water 
surface in each window. This is less likely to be necessary with two waves than with 
one, as the number of unknown parameters is much larger. It should be kept in 
mind, however, as it would be a factor for low order solutions with a large number of 


measurement locations. 
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5.3 Computation Methods 


The LFI method for an array of water surface measurements can be broken down 


into a similar sequence of steps as with the analysis of a pressure record: 
1. Pre-processing of record. 


(a) Determine estimate for level of noise in the record. 
(b) Determine estimate for magnitude and direction of Eulerian current. 


(c) Define a set of continuous records from each gauge from the discrete ob- 


servations by cubic spline interpolation. 
(d) Specify spacing of output locations. 
(e) Compute mean zero crossing frequency. 


(f) Non-dimensionalize record and all parameters. 
2. Primary values of numerical solution parameters are chosen. 


(a) Window width (70) 
(b) Order of solution (J) 
(c) Number of time samples of the water surface records within each window 


(M) 


3. Global solution is computed on an entire wave to provide first estimates for 


local optimization 


4. For each selected output location, a window in the record is defined, and the 
an LFI solution is computed. 
(a) Initial guess for the optimization is determined from the global solution. 


(b) Full nonlinear optimization for all unknown components of the potential 


function is computed. 


(c) Results are checked for spurious solution. 
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i. If no solution, or a spurious solution, is found the solution parameters 


are adjusted, and the optimization repeated. 


ii. If a good solution is found, progress to the next window. 


5.3.1 Pre-Processing of Record 
Accommodating Measurement Error 


Measurement error can be a major source of difficulty with any high order data 
interpretation method. Local methods can be especially sensitive, as each window 
solution relies on the detail contained in a small segment of the record. In applying 
the LFI method to a single pressure trace, Sobey (1992) found it necessary to apply 
a simple moving average filter to field and laboratory data. In his work, the primitive 
kinematic free surface boundary condition was used, requiring an estimate for the 
local gradients in the water surface. In the current work, a modified version of the 
boundary condition (Eq. 5.3) is used which does not require these gradients. This 
makes the method less sensitive to noise, so it was not necessary to apply any filtering 
for the results presented here. 

If there is substantial noise in the measured record, the system of equations can be 
- substantially overspecified, allowing the least squares optimization to accommodate 
the noise in the record. When this is possible, it is preferable to applying a smoothing 
filter to the record, as it does not impose any assumptions on the nature of the record. 
However, if the error bands are very large on the data, it may still be necessary to 


apply filtering to the raw measurements. 


The Mean Water Level 


The mean water depth, h, must be specified as part of the potential function. 
As a time series of the water surface is provided, it is a simple matter to compute a 
mean of the measured records. The mean should be taken over a period much longer 
than a typical wave period, but short enough to accommodate changes in the local 


water level due to astronomical and storm tides. In keeping with the local nature 
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of the approach, this is the local mean water level, rather than a global still water 
level, which might be different, and would be less appropriate. As with subsurface 
measurements, the precision of this value is not critical, as h is only used to locate 


the origin of the coordinate system for the potential function. 


The Eulerian Current 


The Eulerian current, U, and U,, completes the description of the propagation 
medium of the waves, and an accurate estimate of its magnitude and direction is 
important. The measured water surface traces provide no information about the 
current field, so the information needs to be provided from other sources. See Section 


3.3.1 for a detailed discussion of this important parameter. 


Spline of the Records 


As with the previous chapter, a set of continuous records is computed by cubic 
spline interpolation among the measured points at each gauge, allowing complete 
flexibility in the choice of window widths and location of samples in time of the 


records. 


Output Locations 


The spacing of the desired output locations must be chosen to determine the 
placement of the windows. Each window is computed independently so there is no 
restriction on the spacing of the output windows. In addition to selecting the output 
spacing in time, in can be useful when using an array of instruments to determine a 


single central location within the array at which to define the solution. 


Non-dimensionalization 


In order to prevent spurious solutions due to the comparisons of errors of different 
dimensional quantities, all parameters and variables are non-dimensionalized before 


computation with scales defined by the same physically identifiable parameters used 
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in Chapter 3: the mass density of water (g), acceleration of gravity (g), and the mean 


zero crossing frequency of the measured record (w-). 


Length scale = g/w? 


Time scale = 1/w, (5.7) 
3 
Mass scale = ee 
Ww. 


5.3.2 Optimization Procedure 


The nonlinear optimization in the LFI method for the analysis of arrays of water 
surface measurements is very similar to that used for the analysis of a subsurface 
pressure trace. In the single wave approach, the solution is somewhat easier. The 
potential function used by the single wave approach is almost the same as that used 
with a point measurement (both Chapter 3 and Sobey (1992)), with the addition 
of a directional component to the wave number. This results in a single additional 
unknown, but the array of measurements provides at least three times as much data, 
with gauges at at least three spatial locations to specify uniquely the directional 
structure of the sea. The optimization tends to converge more rapidly and robustly 

than with a single point measurement. 
When applying the method with the two wave potential function, there are many 
more unknowns, and the optimization once again becomes somewhat tenuous. As 
with the previous chapter, in is essential to identify a good initial estimate to reduce 


the chances that the optimization will converge to a spurious minimum. 


Global Solution 


The strength of local methods is that they seek a single solution for only a short 
segment of a record at a time. The downside is that a single short segment often 
may not contain sufficient information to identify the directional nature of the wave 
field. In the method presented in the previous chapter the initial estimate could be 
computed from the data in the current window. In contrast, when working with 


spatial arrays, it is necessary to examine a larger segment of the record, for example 
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an entire wave from crest to crest or trough to trough, to get a general sense of the 
directional trend of the wave field. This global solution can then be refined to fit a 
small segment very closely. This approach is used to establish the initial estimate for 
the optimization procedure in each window. 

For the initial estimate to the global solution, it is assumed that the water surface 


records can be approximated with linear wave theory: 
n = acos(kztn + kyyn —wtm +a) (5.8) 
for the single wave method, and 


n = ay cos (een. + kyaYn cae Wit iy a1) (5.9) 


+2 COS (Kr22n =F ky 2Yn i Wotm oP a2) 
for the two wave method, where a, = A,w/g, and A, is the amplitude of the linear 


potential function of the nth wave. 


Directional Trend The first step is to determine the directional trend of the wave 
field. This determination is accomplished by examining the gradients of the water 
surface throughout the segment considered. Estimates for the spatial gradients and 
the elevation of the water surface at the center of the array are computed by finite 
difference approximations. The water surface is expanded in a two dimensional Taylor 


series about the center of the array: 


= One One 
UN(25 2) = ih a Bs) ae (@ = Ba) ae (y — Yc) By +... (5.10) 


where xz, and y, are the coordinates of the center of the array. This expansion is 


written for each of the N measured gauges, at each of the M points in time, resulting 
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in a system of linear equations that can be expressed as the following matrix equation: 


Ll B= 2) We— Mn iol 2 3 o o ina 
If fe=2y Ve = Op 2,1 
Net Ne2 - - - Nem 
One. 9Ne.2 OneM | — 
Ox Ox AW RE Ox 
ONne,1 ONc,2 One,M 
oy Cun so ae ay 
LT Be ay We — hy WN 0 = « « ii 
(5.11) 
O 6 fo) ead 
where zx, and y, are the coordinates of the nth gauge in the array, 7m, se and 
Stent are the water surface and gradients of the water surface at the center of the 


array at time, t,,, and Nr, 1s the water surface elevation at gauge n and timet,,. If 
there are three gauges, the gradients are uniquely specified. If there are more, the 
system is solved in the least squared sense. 

The water surface is traveling either toward or away from the direction of the 
water surface slope depending on whether it is moving up or down at the time. The 
direction is thus determined by the spatial gradient of the water surface, and the sign 


of the time gradient, yielding a set complex direction vectors: 


«(Om \ ( ONm Om 
By = sign (Fe) (Fe + i) (5.12) 


A cubic spline of the water surface at the center of the array (n.) would provide 
an explicit piecewise polynomial expression that could be directly differentiated to 
obtain the gradients in time of the water surface. These estimates would be very 
sensitive to measurement errors in the record. To obtain more robust estimates for 
the time gradients, a smoothing cubic smoothing spline (de Boor 1978) is employed 
instead. This algorithm provides a smoothing parameter, p, that can be set at any 
value between 0 and 1, where p = 0 results in the linear regression fit, and p = 1 
results in the “natural” cubic spline. The smoothing parameter may be varied to 
accommodate varying levels of noise in the record. For the examples in this chapter, 


p = 0.9 was found to be satisfactory. 
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Figure 5.2: Direction vectors for a short crested wave 


For the single wave method, it is expected that the individual directions will vary 
around a central dominant direction, and the propagation direction estimate is the 


orientation of the mean of the complex direction vectors: 


6 = angle (3; SD 3. (5.13) 


For the two wave method, it is expected that the propagation directions at each 
point in time will vary around two dominant directions. In the case of a standing 
wave, these two directions are clearly defined. In a more complex sea, it is not so 
straightforward. To separate the two dominant propagation directions, the set of 
direction vectors, Der is divided into two sets. The two dominant directions are 


defined as the angles of the means of the sets of direction vectors within 7 greater 
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than and less than the mean direction (see Fig. 5.2): 


M, 
1 = e2 
6, = angle (a: ) 3. where: 0< angle( D m) <0+7 
(5.14) 


1 = = 
62 = angle (si SS B. where: 6—-7< angle( Dm) <w) 


Frequency Once the dominant directions have been identified, the approximate 
frequency must be determined. This is accomplished by examining the water surface 
at the center of the array, as interpolated by the finite difference approximation 
described above (Eq. 5.11). The method used is identical to that used in the two 
dimensional method, (Ch. 3): 


LCP Reso 
Mam OF 


(5.15) 


Um = 


On 


O2ieian cm 
5, for the 


where —,3~ is computed from the smoothed spline used for to compute 


direction estimate. If the wave field is unimodal it is expected that there will be a 
single dominant local frequency. The calculated frequencies at each time step will be 
similar in magnitude, and the mean over the record is used as the first estimate of the 
frequency for the single wave. In the two wave method, it is assumed that the bimodal 
sea is the result of reflection, so that the frequency of the incident and reflected modes 
should be the same, and the mean frequency is used as the first estimate for both 


waves. 


Wave Numbers Once the frequency is known, the wave numbers are estimated 


from the linear dispersion relation. 
(w — [Kero = gk, tanh k,h, Kom, = IN COSC, Sion SUK, Sia, — (D-LG)) 


where U,, is the component of the Eulerian current in the direction of the nth wave. 


U, 
eae J2 2 1 Ne )\\ gee 
Un = ,/ U2 + U? cos (an (F ) 6, (5.17) 


xz 
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Amplitudes and Phases The amplitudes and phases of a particular record can 
be found by rearranging the equations as a linear least squares problem by separating 


the cosine and sine components as was done in chapter 3: 


Ne = acos(kzx + kyy — wt + a) 
= bcos (krx + kyy — wt) + esin (kx + kyy — wt) 


(5.18) 
= Ve+ 2 
a@ = arctan (—c/b) 
for the single wave method, and 
Ne = 4, C08 (ket + ky iy — wt + a1) + a2 c08 (ke 22 + ky oy — wt + a2) 
= 6 cos(kpit + kyiy — wt) +c sin (kz iz + ky vy — wt) 
+b» cos (kz,2% + ky2y — wt) + co sin (kz22 + ky 2y — wt) 
a= oi+e; (5.19) 
a= \/b+¢ 
Qa, = arctan (—c;/b;) 
aj = arctan (—c2/b2) 


for the two wave method. The system is determined if 7.(t) is defined at at least two 
points in time for the single wave method, and at four points in time for the two wave 


method. The system is solved in the least squared sense in the case of more points. 


Refining the Linear Estimates These procedures result in very rough estimates 
for the parameters of two intersecting linear waves. The estimates are then refined 


by optimizing for a best linear wave theory fit to the record: 


minimize O(X Sa (Rom = df” (OS Bag Mast oye (5.20) 


n=lem=1 


where: 


f" (X53 2n, Yn tm) = acos (ketn + kyYn — Wtm + @) 
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for the single wave method, and: 


Ff" (X35 lny Yn tm) = a1 C08 (ke En + ky Yn — Witm + 1) 


-++ G2 COS (kioan ain ky,2Yn = Oils tr Ca) 


<= (Rete kya, Qi, kz2; ky, 2,41, a2) 


kn = \V koe a [Be 


for the two wave method. The frequencies are determined from the linear dispersion 


relation: 


Wn = Vgkp tanhk,h + k,Un (5.21) 


where U,, is defined by Eq. 5.17 This optimization results in a linear estimate for the 
water surface that fits the measured records most closely in the segment considered. 

This procedure has been performed on a segment of the records large enough to 
resolve the directional structure of the wave field, usually a complete wave from crest 

-to following crest, or trough to following trough. The final step in computing the 
initial estimates for the final window-by-window optimization is to compute a full 
order global solution to this larger segment of the record. The initial parameters 
for this full order global optimization are provided by the computed linear fit to the 
water surface records, with the amplitudes adjusted for the potential function: 

A, = (5.22) 

Wn 

The higher order Fourier amplitudes are all initially set to zero. 

Using these linear wave theory estimates as the first guess, the full order optimiza- 
tion, (Eq. 5.6), is computed to determine the best full order fit to a global segment of 
the record. The parameters of the potential function computed by this optimization 
are then used as the initial estimate for the final optimization in each defined small 


window in time. 
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Phase Shift The phase parameters, a; and a2 are set for the phase of the global 
segment of the record used. The local time in each window is set so that t = 0 at 
the center of the window. The initial estimate for the phase in each window must be 


shifted to accommodate the change in the time reference frame. 
An = —w,At t+ Qn, global (5.23) 


where At is the difference between the time in the two reference frames. 


Nonlinear Optimization in Windows 


Once the global solution is computed, it is used as the first guess for the parameters 
in each window, and these parameters are refined by solving Eq. 5.6 with a standard 
nonlinear optimization routine. For the results given, the Levenberg-Marquardt algo- 
rithm was used as implemented by the MATLAB Optimization Tool Box (Grace 1992). 
If the optimization successfully converges to a minimum, the solution is checked to 
see if it is a clearly spurious solution. Spurious solutions can be identified by the same 


criteria used in chapter 3: 
e Very large or systematically variable errors. 
e First order amplitude smaller than one of the higher order amplitudes. 
e Unrealistically large or small frequency or wave number. 
e Large discontinuity between the windows in the predicted kinematics. 


As in chapter 3, it is unusual for the routine to converge to a spurious solution. It is 
far more common for the routine not to converge at all. 

If no solution or a spurious solution is found, it is necessary to revise the pa- 
rameters of the solution as was discussed in section 3.3.2. These revisions include 
increasing the width of the window, and if that is not successful, decreasing the order 
of the solution. If neither of these adjustments result in an acceptable solution, the 
window is skipped, and future analysis must be interpolated through that point. As 


with the analysis of a point pressure measurement, these adjustments are most likely 
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to be needed in the long, flat troughs of shallow water waves, or near zero crossings 
in the record. The additional data provided by the array of measurements, and the 
fact that the elevation of the water surface is measured, provides for a much more ro- 
bust optimization than with the subsurface pressure measurement. As a result, these 


adjustments are necessary less often with an array of water surface measurements. 


Locating the Water Surface at the Array Center 


Once the solution is found, the potential function, and thus the complete kine- 
matics in the immediate neighborhood of the array, are defined. The solution is likely 
to be most accurate at the center of the array, and it is often convenient to have a 
solution at a single point, so the water surface at the center of the array must be 
found. This is accomplished by setting up another optimization problem with the 
elevation of the water surface at a few nodes in time throughout the window as the 
unknowns. 

The free surface boundary conditions (Eqs. 5.4 and 5.3) are applied at the hori- 
zontal location of the center of the array, at M points in time throughout the solution 
window. At each point in time, the only unknown is the elevation of the water sur- 
face, and the two boundary conditions provide two independent equations. The water 
surface is defined as that location that results in the least error in the FSBCs, in the 
least squared sense. Each point in time is independent, but the system can be set up 
to solve for a number of points at once. Enough points should be found to specify 
the shape of the water surface throughout the window. For the results given here, 
six points were used in each window (see Figs. 5.6 through 5.23). The mean of the 
measured water surface at the nodes provides a good first estimate. This system 


consistently and rapidly converges to a solution. 
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Figure 5.3: Layout of water surface array 


5.4 Theoretical Records 


5.4.1 Single Wave Records 


The following results are from “measured” data generated by Fourier steady wave 
theory (Sobey 1989). Fourier theory provides a near-exact solution for steady waves 
that provides the complete kinematics. This approach provides a complete set of data 
to compare with the results of the LFI method, without the complications introduced 
by the inevitable errors of data collected in the field or the laboratory. Theoretical 
records were used also because field or laboratory data that included a full set of 
measured kinematics to compare with the results are not available. 

A triangular array of water surface measurements was used, as indicated in Fig. 
5.3. The array is an equilateral triangle with the same dimensions as the DWG-1 
pressure array (Howell 1992). Three measurements were chosen, as three is the min- 
imum number required to provide directional information. Additional measurements 
would provide overspecification, and can easily be accommodated in the formulation. 
With actual field data, additional measurements are recommended, as increasing the 
number of instruments would provide redundancy in the case of instrument failure, 


as well as helping to accommodate measurement error. 
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Kinematics at the predicted water surface 


—— Fourier 
LFl 


Figure 5.4: LFI predictions (J = 3) and exact kinematics at the center of the array 
at the predicted water surface for a steady deep water wave 
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Fig. 5.4 is a steady deep water wave generated by 10th order Fourier theory with 
the following parameters: wave height = 10 m, period = 10 seconds, water depth = 
100 m, zero Eulerian current, and direction of travel 10 degrees from the x-axis. This 
is a fairly large wave in deep water. The window width is 1 sec, 1/10 the zero crossing 
period, and the LFI solution is computed to third order (7 = 3). The LFI method has 
captured the location of the water surface and the kinematics at the surface essentially 
exactly. While linear wave theory might do an adequate job of approximating much 
of the kinematics of a deep water steady wave like this, it is important to remember 
that each of the points in Fig. 5.4 was computed from a small segment of the record 
surrounding that point. In this case, the window width is 1/10 of the zero crossing 
period, or ls. The local nature of this method extends its applicability to irregular 
wave records. 

Fig. 5.5 is a steady shallow wave generated by 18th order Fourier theory with the 
following parameters: wave height = 3 m, period = 10 seconds, water depth = 5 
m, zero Eulerian current, and direction of travel 10 degrees from the x-axis. This is 
a fairly extreme wave in shallow water. The window width is 1 sec., 1/10 the zero 
crossing period, and the LFI solution is computed to third order. As with the deep 
water wave, the LFI method has captured the location of the water surface and the 


kinematics at the surface essentially exactly, including the pronounced sharp crest. 


Choice of Order 


In order to determine the order necessary to accurately capture the kinematics 
of measured waves, it is particularly useful to examine a window near the crest of a 
wave. The crest is usually the region that requires the highest order solution. This is 
particularly true for shallow water waves, but higher order wave theory in all depths 
of water indicates that, as the wave height increases, the crest tends to get sharper, 


and the trough flatter. Capturing this sharp crest requires a high order solution. 


Deep Water A window near a crest of the deep water steady wave given in Fig 5.4 
has been computed at orders 1 through 4. The results in that window are given in 


Fig. 5.6 through 5.13. Figures 5.6, 5.8, 5.10, 5.12 are the non-dimensional errors in 
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Kinematics at the predicted water surface 


Figure 5.5: LFI predictions (J = 3) and exact kinematics at the center of the array 
at the predicted water surface for a steady shallow water wave 
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Crest of a deep water wave 
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Figure 5.6: Free surface boundary condition errors for a deep 
water wave at order 1 
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Figure 5.7: Free surface and velocities at the center of the array 
for a deep water wave at order 1 
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Crest of a deep water wave 
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Figure 5.8: Free surface boundary condition errors for a deep 
water wave at order 2 
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Figure 5.9: Free surface and velocities at the center of the array 
for a deep water wave at order 2 
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Crest of a deep water wave 
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Figure 5.10: Free surface boundary condition errors for a deep 
water wave at order 3 
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Figure 5.11: Free surface and velocities at the center of the array 
for a deep water wave at order 3 
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Figure 5.12: Free surface boundary condition errors for a deep 
water wave at order 4 


Crest of a deep water wave 


04 -03 -0.2 -0.1 0 0.1 0.2 0.3 0.4 


0.1 
— actual 
WW, 6 © predicted 
g 


04 -03 -02 -0.1 (0) 0.1 0.2 0.3 0.4 


Figure 5.13: Free surface and velocities at the center of the array 
for a deep water wave at order 4 
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the free surface boundary conditions (Eqs. 5.3 and 5.4) at each of the measured nodes. 
These are the data that would be analyzed in a practical situation. Figures 5.7, 5.9, 
5.11, 5.13 are a comparison between the predicted and actual values for the water 
surface and the velocities at the surface at the center of the array. The actual values 
were computed using Fourier wave theory. The actual values would not be available 
for comparison when analyzing field records. 

The first order computation results in errors in the free surface boundary condi- 
tions of only order 10~° (Fig. 5.6) as well as very accurately predicting the velocities 
at the water surface at the center of the array (Fig. 5.7). It is a surprisingly accu- 
rate first order solution. This is because of the local nature of the method. When 
a single first order solution is used to capture the entire wave, the error is larger, of 
order 107°. It also should be noted that this first order solution is not the same as a 
linear wave theory solution, even locally. The full nonlinear boundary conditions are 
preserved, and the frequency and wave number are free to vary, and are not bound 
by a dispersion relationship. For steady waves in deep water, linear wave theory is 
fairly accurate. Linear theory is not, however, a local solution, and is not directly 
applicable to irregular waves. 

At second order, the free surface boundary condition errors are smaller, of order 
10-®, and the velocities at the surface match the Fourier solution visually perfectly. 
At third and fourth order, the errors in the free surface boundary conditions continue 
to decline, and the water surface velocities continue to match the Fourier solution 
well. In deep water, for waves of this height, second order is more than adequate to 
capture the surface kinematics of this wave. Higher order solutions are likely to be 
necessary to capture the irregularity of field records, even in deep water. 

Choice of order is dictated by the desired accuracy of the solution, and by the ease 
of convergence to the solution. As there are more free parameters at higher order, 
a higher order solution will always have smaller errors in the free surface boundary 
conditions. In the case of this example, the solution converged at all orders very 
quickly, so there is little penalty in using third or fourth order. With an irregular 
record, in contrast, convergence can be more difficult, and it is occasionally necessary 


to resort to lower order. 
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As the order is increased, there are more free parameters, and care must be taken 
to include sufficient points in each window. As discussed above, for fairly low orders 
and an array of measurements, the limiting factor is the minimum number of points 
needed to define the curvature in the window, rather than to specify the system of 
equations. The above examples were computed from an array of three points, and 
so required a minimum of only one sample to specify the system at first and second 
order, and two points for up to order 8. When attempted, the solution would not 
converge with only one sample. With two samples, the optimization algorithm found 
a reasonable solution, but this small number of samples would not define the shape of 
the water surface adequately if it were part of an irregular record, and so three points 
were used in each window for all orders. In general, a minimum of three points should 
be used, and more may be necessary to accommodate a highly irregular profile, or 
the inevitable measurement error in a field record. 

In the case of theoretically generated records the need for more sampling of points 
poses no problem, but with field records, there are limitations as to the spacing of 
the sampled points. In order to free up the spacing of points for the LFI method, 
points are sampled from a cubic spline interpolation of the actual record. This allows 
the points to be sampled anywhere within the window. While computationally it is 
' possible to sample as many points as necessary in a small window, if that window is, 
in fact, defined by only a couple of actual data points, it is not appropriate to try 
to fit a high order solution to a segment defined by only a few observational points. 
In order to include sufficient actual data points to justify the increased order, the 
window must be increased in size. While increasing the size of the window permits a 
higher order solution, it also compromises the local nature of the method. The goal 
of the LFI method is for the solution to be as local as possible, which is achieved by 


selecting as small a window as possible at fairly low order. 


Shallow Water A window near a crest of the shallow water steady wave given 
in 5.5 has been computed at orders 1 through 5. The results are given in Fig. 5.14 
through 5.23. These figures are analogous to those previously discussed, but on a 


shallow water wave. 
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Figure 5.14: Free surface boundary condition errors for a shallow 
water wave at order 1 
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Figure 5.15: Free surface and velocities at the center of the array 
for a shallow water wave at order 1 
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Figure 5.16: Free surface boundary condition errors for a shallow 
water wave at order 2 


Crest of a shallow water wave 


— actual 


-0. 
02 © predicted 
WW» 
0 4 
eae) 


SAT Ae -0.1 (0) 0.1 0.2 0.3 0.4 


Figure 5.17: Free surface and velocities at the center of the array 
for a shallow water wave at order 2 
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Figure 5.18: Free surface boundary condition errors for a shallow 
water wave at order 3 
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Figure 5.19: Free surface and velocities at the center of the array 
for a shallow water wave at order 3 
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Figure 5.20: Free surface boundary condition errors for a shallow 
water wave at order 4 
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Figure 5.21: Free surface and velocities at the center of the array 
for a shallow water wave at order 4 
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Figure 5.22: Free surface boundary condition errors for a shallow 
water wave at order 5 
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Figure 5.23: Free surface and velocities at the center of the array 
for a shallow water wave at order 5 
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In this case, the first order computation results in substantial errors in the free 
surface boundary conditions of order 10~? (Fig. 5.14) and does not predict the ve- 
locities at the water surface at the center of the array as well as it did in deep water 
(Fig. 5.15), noticeably under-predicting the horizontal velocities. The local nature 
of the method allows a first order solution to get fairly close, but steady waves in 
shallow water are highly nonlinear, and a first order solution simply can not capture 
the sharp crest. 

At second order, the free surface boundary condition errors are slightly smaller, 
but still of order 10~?, and still under predicting the horizontal velocities at the 
surface. At third and fourth order, the errors in the free surface boundary conditions 
continue to decline, but the water surface velocities do not match the Fourier solution 
visually perfectly until fifth order (Fig. 5.23). In shallow water, it may be necessary 
to use up to fifth order to capture accurately the surface kinematics of a fairly large 
wave. In the case of this example, the solution converged at all orders very quickly, 
and there is little penalty in using fourth or fifth order. With an irregular record, 
convergence will be more difficult, and it may sometimes be necessary to resort to 
lower order. 

As previously discussed, the use of higher order solutions requires that more points 
’ be sampled in each window. This is not likely to be a limitation with the single wave 
form with an array of measurements, as it is likely to be necessary to overspecify the 


system in order to define the curvature within the window. 


5.4.2 Records of Two Intersecting Waves 


In order to develop the method using two intersecting waves, theoretical water 
surface records were generated by Ohyama’s fourth order intersecting wave theory 
(Ohyama, Jeng, and Hsu 1995a). Ohyama’s method is a Stokes-style solution for 
irrotational intersecting waves that is accurate to fourth order, and applicable in deep 
water (see Appendix A). As the resulting water surface from intersecting waves can 
be quite complex, a desired wave height can not be specified. Rather, the amplitude 


of the first order component of each of the intersecting free waves is specified. The 
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first order components are the largest components of the resulting wave field, and thus 
give an approximation of the size of the resulting waves. This wave theory assumes 


a zero Eulerian current. 


Standing Wave Figure 5.24 shows the results of the method for a standing wave 
generated by two identical intersecting waves traveling in opposing directions. The 
parameters of the waves are: Period = 10 sec., first order amplitudes = 3 m, prop- 
agation directions, 15 and 195 degrees from the x axis, in deep water. The LFI 
method to third order (J = 3) finds the kinematics at the water surface almost ex- 
actly. While these results show the complete wave, it is important to keep in mind 
that as with all the previous results, each of the indicated points is in the center of 
a separate window, and was computed independently of the other windows. In this 
case, the standard window width was 1 sec., or one tenth the zero-crossing period of 
the wave, with a third order potential function (J = 3), and four water surface nodes 
distributed equidistantly in time in each window (M = 4), at each gauge, for a slight 


overspecification, with 24 equations in 20 unknowns 


Short Crested Wave Figure 5.25 is the water surface at a point in time of the 
short crested wave that would result in the reflection of a steady wave from a sea 
wall, as indicated by figure 5.26. The parameters of the wave are: Period = 10 sec., 
first order amplitudes = 3 m, propagation directions, 30 and 150 degrees from the x 
axis, in deep water. Figure 5.27 shows the results of the method for this short crested 
wave. The LFI method to third order again finds the kinematics for this wave almost 
exactly. In this case, the standard window width was also 1 sec., or one tenth the 
zero-crossing period of the wave, with a third order potential function (J = 3) and 
four water surface nodes (M = 4) in each window, for a slight overspecification, with 


24 equations in 20 unknowns. 


Choice of Order 


As with the previous examples with steady waves, it is useful to examine a window 


near the crest in order to determine the order necessary to accurately capture the 
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Figure 5.24: LFI predictions (J = 3) and exact kinematics at the center of the array 


at the predicted water surface for a standing wave 
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Figure 5.26: Schematic of reflection that would generate the short crested wave given 


in figure 5.25 


kinematics of a short crested wave. 

A window near a crest of the short crested wave given in Fig. 5.27 has been com- 
puted at orders 1 through 4. The results in that window are given in Fig. 5.28 through 
5.35. Figures 5.28, 5.30, 5.32, 5.34 show the errors in the free surface boundary con- 
ditions at each of the measured nodes. These are the data that would be analyzed in 
a practical situation. Figs. 5.29, 5.31, 5.33, 5.35 give a comparison between the pre- 


dicted and actual values for the water surface and the velocities at the surface at the 
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Figure 5.27: LFI predictions (J = 3) and analytical kinematics at the predicted water 


surface for a short crested wave 
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Figure 5.28: Free surface boundary condition errors for a short 
crested wave at order 1 
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Figure 5.29: Free surface and velocities at the center of the array 
for a short crested wave at order 1 
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Figure 5.30: Free surface boundary condition errors for a short 


crested wave at order 2 
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Figure 5.31: Free surface and velocities at the center of the array 


for a short crested wave at order 2 
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Figure 5.32: Free surface boundary condition errors for a short 
crested wave at order 3 
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Figure 5.33: Free surface and velocities at the center of the array 
for a short crested wave at order 3 
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Figure 5.34: Free surface boundary condition errors for a short 
crested wave at order 4 
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Figure 5.35: Free surface and velocities at the center of the array 
for a short crested wave at order 4 
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center of the array. The theoretical values used Ohyama’s fourth order intersecting 
wave theory. The actual values would not be available for comparison when analyzing 
field records. 

It is clear that the first order computation results in substantial errors in the 
free surface boundary conditions of order 107° (Fig. 5.28), as well as significantly 
underestimating the horizontal velocities (Fig. 5.29). As with the steady deep water 
wave previously discussed, the first order solution is quite reasonable. This is because 
of the local nature of the method, and the fact that it is not a linear solution. The full 
nonlinear boundary conditions are imposed, and the frequencies and wave numbers 
are free to vary, and are not bound by a dispersion relationship. 

At second order, the free surface boundary condition errors are better, of order 
10-°, and the velocities at the surface match the analytical solution visually perfectly. 
At third and fourth order, the errors in the free surface boundary conditions continue 
to decline, and the water surface velocities continue to match the theoretical solution 
well. 

As with the single wave method, choice of order is dictated by the desired ac- 
curacy of the solution, and by the ease of convergence of the optimization. For the 
above examples, four water surface nodes (M = 4) distributed equidistantly in time 
in each window, at each gauge, at orders 1 through 3. These values providing an over- 
specification and sufficient points to define the shape of the water surface throughout 
the window. Five points (M = 5) were used at fourth order, as four points provide 
only 24 equations, and there are 28 parameters to be found. Five points provides 30 
equations for a slight over specification. At fifth order, there are 38 unknowns, and 
seven points in time (M = 7) would have to be used in each window. Sampling this 
many points in a single window would require very closely spaced data or a wider 
window. Another solution would be to use an array with additional measurement 
locations. An array of four gauges, for example, would provide eight equations per 
point in time, and would allow a fifth order solution to be computed from five points 
(M = 5) per window. It would be advantageous to use as many gauges as possible in 
shallow water, where higher order solutions are necessary. 


Unfortunately, the theoretical solution used for this analysis, while the best avail- 
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able for intersecting waves, is only accurate to fourth order, and only being applicable 
in deep or transitional water. This precluded an analysis of the behavior of the so- 
lution in shallow water. When the method is applied in shallow water, the lessons 
learned in the two dimensional case should be applied, and a higher order solution, 
perhaps to fifth order, should be used. 

It is also the case that the Ohyama solution is a single solution that is accurate 
to fourth order globally. The local nature of the LFI method allows each separate 
window in time to be represented by a unique solution, representing the entire record 


with far more accuracy than a global solution of the same order. 


Choice of Window Width 


As mentioned in the previous section, it is important to keep the window width at 
a minimum. For the short crested wave shown in Figure 5.27, the standard window 
width was one tenth the zero-crossing period (0.17,). It was necessary to widen the 
window to 0.157. at three locations in order to find a solution. This indicates that 
the standard window width is set close to the smallest size that would be effective. 
In general, it is necessary for the window to contain enough of the record to define 
the local curvature, and include sufficient data points to reasonably expect to define 
a high order solution. Most often, the limiting factor will be the sample rate of the 
data. For example, the above examples are computed for a short crested wave with 
a period of ten seconds. In those examples, the window width of 1/10 the period 
(1s) was adequate to capture most of the wave. Unfortunately, field records are often 
sampled at only 1Hz. At this sampling frequency, and a one second window, there 
would be a maximum of only two data points in each window, and it may be necessary 


to widen the window to include more data. 


Effects of Array Size 


The examples above are all computed from data sampled by a fairly small array 
(Fig. 5.3). The 1.6m size of the array is about 0.01 of the wavelength of 10s wave in 


deep water. That particular size was chosen because it has the same dimensions as 
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the DWG-1 pressure array that is used for the analysis of pressure records in Ch. 6. 
A small array has the advantage of having a small distance to interpolate to the 
center of the array, and is local in the sense that each of the sensors is in essentially 
the same part of the wave. The disadvantage of a small array is that the differences 
between the measured values at each of the nodes are small, and thus more sensitive 
to measurement error. With the theoretical records used for the previous analysis, 
the method was not sensitive to array size, giving equally accurate results with arrays 
up to ten times the size of the DWG-1. The spacing between the gauges of an array 
will determine the wavelengths to which the gauge is most sensitive. Window widths 
of approximately one tenth a typical period have been found to be effective. This 
indicates that a gauge spacing of about one tenth of the expected dominant wave 


length might be an optimum spacing. 


5.5 Laboratory Measurements 


The results from the analysis of analytically derived records given in the previous 
sections demonstrate the method’s capabilities in a variety of conditions. They have 
also been very useful in helping to determine the range of solution parameters that 
must be adopted, including order of solution, window width, and the number of time 
samples taken in each window. The question still remains, however, as to how well 


the method works with actual irregular waves. 


5.5.1 Laboratory Experiment 


The following results are taken from a laboratory experiment performed together 
with Dr. Steven A. Hughes at the Coastal Engineering Research Center, Army Corps 
of Engineers Waterways Experiment Station in Vicksburg, MS, in June of 1995. The 
waves were generated in a 0.46m wide flume by a programmable, hydraulically driven, 
piston type wave board (see Fig. 5.36). At the other end of the flume was a beach 
with 1:30 slope, and a single layer of wave-absorbing horsehair matting. This beach 


has been shown to have bulk reflection coefficients that vary from about 5% for 
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Figure 5.36: Schematic of wave flume 


Hp»0/h = 0.1 up to about 15% for Hmo/h = 0.4 (where Hmo is the incident spectral 
significant wave height and h is the water depth) (Sultan 1992; Sultan and Hughes 
1993). 

The water surface was measured by capacitance wave rods, calibrated with a cubic 
calibration function. The velocity data were collected using a Dantek laser Doppler 
velocimeter (LDV) system operated in the backscatter mode. The LDV system fea- 
tured a 2-watt argon-ion laser equipped with a fiber-optic probe that measures two 
orthogonal water velocity components (horizontal and vertical). Velocity data were 
converted in real time to engineering units (m/s) and written to a computer file 
simultaneously with the wave rod data. 

The wave rods and LDV were placed near the middle of the flume, with the wave 
rods arranged in an equilateral triangle with leg length of 0.17m (see Fig. 5.37). The 
LDV was situated to measure the horizontal and vertical velocities at the center of the 
array, at a variety of vertical elevations. The flume was allowed to reach quiescence 


in between runs, and the waves were measured for only a short time after starting the 


118 


e Wave maker ———~> 


— 


LDV location 0.17m 


gauge 3 gauge 1 
Plan 


Figure 5.37: Layout of array in flume 


wave maker. At a depth of 0.35m, the approximate wave celerity (\/gh) is 1.8m/s. As 
the measuring location is located 9m from the beach, it takes about 10s for the waves 
to travel from the measurement location to the beach, and for the reflected waves 
to travel back to the study location. Thus, the first 10s of the record are assured of 
being uncontaminated by reflection from the beach. Using only the beginning of the 
record also prevents a Stokes drift induced return current to be established, and so 


the Eulerian current is taken to be zero. 


5.5.2 Laboratory Results 


Figure 5.38 shows a sample wave taken from an irregular wave record measured 
in the flume. The waves were generated with a target TMA spectrum (Bouws et al. 
1985) with a peak period of 2.5s and mean wave height of approximately 12cm. The 
water depth was 0.35m, with the velocities measured by the LDV at an elevation of 
0.05m below the mean water level. The solid lines on figure 5.38 are the horizontal and 
vertical velocities as measured by the LDV. The circles are the velocities predicted by 
the LFI method at the center of each computed window. The solid line in the water 
surface plot is the measured elevation at wave rod 2, at the same x coordinate as the 


center of the array. The circles are the water surface elevations at the center of the 


Figure 5.38: 


Water surface 


Results on a sample wave in a flume 
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Figure 5.39: Free surface boundary condition errors at the crest 
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Figure 5.40: Computed and measured velocities near the crest 


of a sample wave in a flume 
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array as computed by the LFI method. 

This solution was computed to third order, with a window width of 1/10 the mean 
zero crossing period of the record (approximately 0.2s). Four time samples were taken 
in each window, providing substantial overspecification of the system, and allowing 
the shape of the water surface to be well defined in each window. The first trough is 
fairly flat, and the optimization did not converge to a solution quickly with the small 
window. After the window was widened, the optimization converged quickly to the 
given results. 

The LFI method has captured the detail of the kinematics of this irregular wave 
very well. The comparisons to the measured data are given just below the troughs of 
the wave. Note that the LFI method accurately captured the secondary hump in the 
vertical velocity (w) near the 2.5s point. The data used to compute the kinematics 
were the measured water surface. The computed water surface was found from the 
predicted kinematics, and it matches the measured water surface almost exactly. The 
accurate computation of the water surface indicates that the predicted kinematics 
near the surface may, in fact, be more accurate than the predictions at the elevation 
where the kinematics were measured. This is to be expected, as the data used to 
compute the solution was measured at the surface, and thus the results should be 
- most accurate there. Achieving the greatest accuracy near the surface can be useful, 
as the surface is the location where the water velocities and accelerations are greatest. 

Figures 5.39 and 5.40 show the results of the computation in a window near the 
crest of the wave given in Fig. 5.38. The errors in the DFSBC are of order 107%, and 
substantially smaller in the MKFSBC. The computed and measured velocities do not 
match exactly, but the magnitude and trend are very close. The results are similar 


in the other windows. 


5.6 Discussion 


The LFI method is an effective and efficient way to re-create the kinematics of 
irregular waves from the measurements of an array of wave rods. It is able to capture 


the kinematics from primarily uni-directional seas as well as the seas near a reflect- 
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ing surface. The local nature of the approach allows a nonlinear solution without 
prohibitive computational costs, using a fairly simple form for the potential function, 
and allowing the parameters of the potential function to change with time. 

The examples in this Chapter were all computed using MATLAB running under 
the Linux operating system on an Intel 9)MHz Pentium processor based PC. The 
shallow water steady wave (Fig 5.5) took the longest to compute, about 30min, and 
85x 10° floating point operations (flops). The steady deep water wave (Fig. 5.4) took 
the least time, 10min. and 15x10® flops, the standing wave (Fig. 5.24) took 20min. 
and 244x10® flops, the short crested wave (Fig. 5.27) took 11min. and 93x10® flops, 
and the laboratory record (Fig. 5.38) took 21min. and 83x10® flops. 

The examples given in this chapter provide guidance as to the parameters of a 
solution to be applied to field records. Far from reflecting surfaces, using a potential 
function representing a single wave is effective. Near a reflecting surface, a potential 
function representing two nonlinear intersecting waves is capable of capturing the 
standing or short crested waves that are likely to develop. 

Window widths of 1/10 of the zero crossing period are small enough to maintain 
the local nature of the solution, and capture the detail of the record, while being 
large enough to include the local trend of the wave field. On certain segments of the 
record, the window width must be increased in order to include sufficient curvature 
in the record to find a solution. The widening of the window is most often needed 
in long, flat troughs in shallow water, or near zero crossings of the record. In either 
case, the window rarely needs to be larger than 1/5 of the zero crossing period. 

Low order solutions are quite adequate for deep water waves. For most waves, 
second order is adequate. For the very steepest waves, slightly higher order may be 
appropriate, and there is little computational penalty in including the higher order 
terms. In shallow water, higher order solutions are necessary. Third order is adequate 
in many cases, but including up to fifth order is recommended for extreme waves. 
When attempting a high order solution such as this with the two wave method, it 
may be necessary to use arrays with more than three gauges in order to provide 


sufficient equations to specify the solution. 
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Chapter 6 


Array of Subsurface Pressure 


Measurements 


Arrays of subsurface pressure gauges are installed in the sea in order to capture 
directional wave field data. These instruments have an advantage over arrays of 
water surface wave rods, as they are mounted below the surface, where they are less 
susceptible to damage, by either the strong wave forces experienced during storms, 
as well as vandalism or accidents with vessels. 

Unfortunately, the reason they are less susceptible to damage from wave forces is 
because they are placed under the surface, often at the sea bed, where the action of 
waves is the smallest. This leads to the mathematically ill-posed problem discussed 
in chapter 3. This difficulty is somewhat mitigated by the added information made 
available by the multiple gauges in the array, but still limits the ability to determine 
the detail of the kinematics, particularly when there is a lot of high frequency energy 


in the sea state. 


6.1 Formulation of Solution 


The lessons learned in the previous chapters provide for a fairly straightforward 
determination of the formulation for adapting the LFI method to the analysis of 


arrays of pressure gauges. The formulation for a single subsurface pressure gauge was 
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presented in chapter 3, and for arrays of water surface measurements in chapter 5. 
The formulation for arrays of subsurface pressure gauges is very similar to that for 
arrays of water surface measurements, with the addition of the need to determine the 
location of the water surface, as was done with the analysis of a subsurface pressure 
trace. The following description includes the required information; greater detail has 
been given in previous chapters. 

The flow is assumed to be irrotational and incompressible, with a potential func- 
tion that represents either a single directional wave, or two separate intersecting 
waves. The potential function locally representing a single wave is Eq. 4.10 reduced 
to a single wave: 
cosh 7A (h + z) 

cosh j Ah 


k= eae 


where U, and U, are the components of the known depth uniform Eulerian current 


i 
(x, z,t) =Urzr + Uyy + > A; 


j=l 


sinj(kzr +kyy—wt+a) (6.1) 


in the z and y directions, h is the mean water depth, J is the truncation order of the 
Fourier series, A; are the Fourier coefficients, w is the local fundamental frequency, 
k, and k, are the components of the local fundamental wave number in the z and y 
directions, and A is the magnitude of the local wave number. 

The potential function for two intersecting waves to fourth order is Eq. 4.18, 


repeated here: 


20 
oz, Y, 2; t) = U;,z + Uyy ale S- A; C(K;h, K;z) sin 0; (6.2) 
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At first order: 


a1 = (51) 
At second order: 
G3 = (2Sa)) 
05 = (5; + S2) 
At third order: 
07 = (35) 


09 = (2S, + S2) 
oy, = (Si + 282) 


And at fourth order: 


o13 = (451) 
o15 = (35; + S2) 
O17 = (25; + 25>) 
o19 = (51 + 352) 
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oD = (S2) 
04 = (252) 
06 = (Sy — 52) 
0g = (352) 


010 = (2S, = 52) 
012 = (S) = 2S) 


Both of these potential functions solve the field equation (Laplace, Eq. 4.2) and the 
bottom boundary condition (Eq. 4.3) exactly. The rest of the solution is determined 


by the free surface boundary conditions: the modified free surface boundary condition 
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the dynamic free surface boundary condition (f?), 
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In order to apply the free surface boundary conditions, the location of the free 
surface must be found, together with the potential function, in each window. In 
order to locate the free surface, the water surface is defined at M nodes, located at 
the horizontal location of the center of the array, and equidistantly spaced in time 
throughout the window. The elevation of these nodes is unknown, and will be sought 


as part of the solution. 


6.2 Formulation of the Optimization 


The formulation of the optimization for the LFI method as applied to the inter- 
pretation of an array of pressure measurements has a great deal in common with the 
formulation for an array of water surface measurements and the formulation for a sub- 
surface pressure record (chapters 3 and 5). Only the framework is presented, together 


with any additional information unique to the application to a pressure array. 


Observational Equations The observational equations are the equations that con- 
fine the solution to fit the measured records. The predicted subsurface pressure is 
available from the Bernoulli equation, Eq. 6.5, which is applied at the location of the 
gauges in the array to define the observational equations. The error in the Bernoulli 
equation is the difference between the measured dynamic pressure at the given lo- 


cation and time and the dynamic pressure predicted from the kinematics defined by 
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the potential function. Sufficient independent equations are defined by applying the 
Bernoulli equation at a number of points in time throughout the considered window. 

In order to specify the solution to a system of equations, there must be at least 
as many independent equations as unknown parameters in the system. In order 
to specify the solution to the LFI formulation for a pressure array, the boundary 
conditions (f? and f*) are applied at each of the water surface nodes, yielding 2M 


B.) is applied at J times on each of the N 


n,t 


equations, and the Bernoulli equation ( 
measured pressure records within the local window, yielding NJ equations. In the 
single wave case, there are 4+ J + M unknown parameters sought in Eq. 6.1 (kz, ky, 
w, a, Ay... Aj, and 7...) in each window, so that if (2M@+N/)>(4+J+M) 
the solution is specified. In the two wave case, there are 2)> 7 +8+M unknowns in 
Eq. 6.2 (2k,, 2ky, 2w, 2a, Ay... Az, and m...nvw: 10+ M at first order, 14+ M at 
second order, 20+ M at third order, and 28+ M at fourth order). The solution is 
specified when (2M + NI) > (25>7 +8+M). This system results in the following 


nonlinear least squares optimization. 


minimize O(X yr DO aa gin ae 
n=1 t=1 (6.7) 


SPX: Ley Ye, "Mm; aa) ap fn (Xs Le, Ye, Mm, bia) 


m=1 


where: 
28 = (Re, py Oy Oh Bigg coe g Agi o0- Mn) 
for the one wave case, and: 
DC (eet pil she Cig [eobin | Way Ob, As oo guava Mil os olfhve) 


for the two wave case. ,, Yn, and z,,p are the coordinates of the nth gauge, x. and y. 
are the coordinates of the center of the array, and 77,, is the elevation of water surface 
node m. 

As with the previous analysis, overspecification is helpful in accommodating the 
measurement errors in field records, as well as being required to define the shape of 


measured records and water surface in each window. 
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6.3 Computation Methods 


The LFI method for an array of pressure measurements can be broken down into 


the same sequence of steps as with the analysis of a water surface array: 
1. Pre-processing of record. 


a) Determine estimate for level of noise in the record. 


b 


Determine MWL and subtract hydrostatic pressure from the records. 
d 


(a) 
(b) 
(c) Determine estimate for magnitude and direction of Eulerian current. 
(d) 


Define a set of continuous records from each gauge from the discrete ob- 


servations by cubic spline interpolation. 
(e) Specify spacing of output locations. 
(f) Compute mean zero crossing frequency. 


(g) Non-dimensionalize record and all parameters. 
2. Primary values of numerical solution parameters are chosen. 


(a) Window width (79) 

(b) Order of solution (J) 

(c) Number of time samples of the pressure records within each window (1) 
(d) Number of water surface nodes within each window (M) 

3. Global solution is computed on an entire wave to provide first estimates for 


local optimization 


4. For each selected output location, a window in the record is defined, and an 
LFI solution is computed. 
(a) Initial guess for the optimization determined from the global solution. 


(b) Full nonlinear optimization for all unknown components of the potential 


function is computed. 


(c) Results are checked for spurious solution. 


i. If no solution, or a spurious solution, is found, the solution parameters 


are adjusted, and the optimization repeated. 


ii. If a good solution is found, progress to the next window. 


6.3.1 Pre-Processing of Record 


The pre-processing of the measured records is essentially the same as presented in 
the previous chapters. Estimates for the mean water level (h) and Eulerian current 
(U, and U,) must be computed to define the propagation medium. Cubic splines of 
the measured records are computed to provide continuous records for computation, 
and the desired output locations are chosen. The records and all parameters and 
variables are non-dimensionalized by the following parameters: the mass density of 
water (p), acceleration of gravity (g), and the mean zero crossing frequency of the 


measured records (w,). 


6.3.2 Optimization Procedure 


The nonlinear optimization in the LFI method for the analysis of arrays of pressure 
measurements is very similar to that used for the analysis of a single subsurface 
pressure trace. In the single wave approach, the solution is somewhat easier. The 
potential function used by the single wave approach is almost the same as that used 
with a point measurement (both chapter 3 and Sobey (1992)), with the addition 
of a directional component to the wave number. This results in a single additional 
unknown, but the array of measurements provides at least three times as much data, 
with gauges at at least three spatial locations to specify uniquely the directional 
structure of the sea. The result is that the optimization tends to converge more 
rapidly and robustly than with a single point measurement. 

When applying the method with the two wave potential function, there are far 
more unknowns, and the optimization once again becomes somewhat tenuous. As 


with the previous chapters, it is essential to identify a good initial estimate to reduce 
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the chances that the optimization will converge to a spurious minimum. 


Global Solution 


As with arrays of water surface measurements, a single short segment of the mea- 
sured pressure records often may not contain sufficient information to identify the 
directional trend of the wave field. As a result, it is necessary to examine a larger 
segment of the record, for example an entire wave from crest to crest or trough to 
trough, to get a general sense of the directional trend of the wave field. This global 
solution can then be refined to fit a small segment very closely. This approach is used 
to establish the initial estimate for the optimization procedure in each window. 

The procedure for computing the initial estimate to the global solution for an 
array of pressure measurements is virtually identical to that used in the previous 
chapter for arrays of water surface measurements. It will be outlined here. For more 
detail, see section 5.3.2. 

For the initial estimate, it is assumed that the pressure records can be approxi- 


mated with linear wave theory: 
Pa = acos(kzt, + kyYn — wtm + a) (6.8) 


for the single wave method, and 


Pa = 4 C0S(ke itn + ky Yn — Witm + A) (6.9) 
+a2 C08 (kz,22n + ky,2Yn — Watm + Q2) 
for the two wave method, where: 
an = AnpWn coy Za as *») 
cosh kh (6.10) 


Ky, = \/k?2 + k2 


A, is the amplitude of the linear potential function of the nth wave, and <z, is the 


elevation of the pressure array. 
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Directional Trend A first estimate for the directional trend of the wave field is 
computed from the gradients of the dynamic pressure. A set of complex direction 


vectors are computed: 


D, =sien (“ix ) (“B 4 item.) (6.11) 


The spatial gradients are estimated from finite difference approximation from the 


measured pressures, and the time gradient is computed from a fitted smoothing spline 
(de Boor 1978) of the estimated dynamic pressure at the center of the array. 
For the single wave method, the estimated propagation direction is the orientation 


of the mean of the complex direction vectors: 


1 
6 = angle (3s ZX 3 (6.12) 


For the two wave method, the two dominant propagation directions are defined 
as the angles of the means of the sets of direction vectors within 7 greater than and 


less than the mean direction: 


M 
6, = angle (a SS 3. where: 0< angle( Dm) <0+4+7 
(6.13) 


M> 
6. = angle (x SS 3. where: Oa=mr< angle( Dm) <6 


Frequency The frequency is computed from the second derivative of the pressure 
at the center of the array. 


1 O-pitem 


ee 6.14 
Pd,cjm ot? ( ) 


Wm = 


8? Da.c.m 


at? 
wave method, it is assumed that the bimodal sea is the result of reflection, so that 


where is computed from the same smoothed spline used above. For the two 


the same frequency is used as the first estimate for both waves. 


Wave Numbers The wave numbers are estimated from the linear dispersion rela- 


tion. 


(w — K,U,)? = gK, tanh K,h, [Big cy IS  COS Oa. OG sim ae (Gell) 
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where U,, is the component of the Eulerian current in the direction of the nth wave. 
U, = ,/U2 + U2 cos (can (=) - é) (6.16) 


Amplitudes and Phases The amplitudes and phases the record are found by 
separating the cosine and sine components so that the equations can be solved as a 
linear least squares problem. 
Dac = acos(k,x + kyy — wt + a) 
= bcos (k,z + kyy — wt) + csin (k,x + kyy — wt) 


(6.17) 
G=ViPEe 
a@ = arctan (—c/b) 
for the single wave method, and - 
Dac = 4,008(ke it + kyiy — wt + a1) + a2 008 (ker + ky oy — wt + a2) 
= 6, cos(kzix + kyiy — wt) +c, sin (kz iz + ky iy — wt) 
+-b2 cos (kz,22 + ky,2y — wt) + co sin (kz2r + ky2y — wt) 
Th = b? + c? (6.18) 
a= /+c 
Q, = arctan (—c,/b;) 
Qa, = arctan (—c2/b2) 


Refining the Linear Estimates These rough estimates are refined by optimizing 


for a best linear wave theory fit to the record: 


N M 
minumize O(X) = DS Ni (ives 3 22 OS Monta) (6.19) 


(=) ToSil 


where: 


fP4 (X35 Ln, Yn, tm) = acos (kptn + kyYn — wtm + @) 


X = (kz, ky, a, a) 


K=,/k2+k 
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for the single wave method, and: 


fP4 (X35 Ln, Yns tm) = a1 008 (kein + ky i Yn — Witm + 1) 


+ a2 608 (kz22n + ky,2Yn — Watm + Q2) 


xX = (ke, kya, 1, kz2, ky2, a2, 41, a2) 


Kn = [Rn + hey 


for the two wave method. The frequencies are determined from the linear dispersion 


relation: 


wn = /gK, tanh K,h + K,U, (6.20) 


where U,, is defined by Eq. 6.16 This optimization results in a linear estimate for 
the dynamic pressure that fits the measured records most closely in the segment 
considered. 

The final step in computing the initial estimates for the final window-by-window 
optimization is to compute a full order global solution to this large segment of the 
record. The initial parameters for this full order global optimization are provided by 
the computed linear fit to the water surface records, with the amplitudes adjusted for 
' the potential function: 


An cosh h,,h 
7 poh coshlka(h-= 25) 


(6.21) 


The higher order Fourier amplitudes are all initially set to zero. 


Water Surface The location of the water surface at M nodes in the center of the 
array throughout the segment is estimated from the linear pressure response function 


with stretching (Nielsen 1989): 
a, (tm) cosh (k (h + pa(tm)/pg)) 


Wiphtre = RETR Aaa 


= (6.22) 

pg cosh k(h + zp) 
where k is the mean of the two wave numbers, Zp is the elevation of the pressure 
array, and 7, is the elevation of the water surface node at the center of the array at 


tm- Pa(tm) is computed from Eq. 6.8 or Eq. 6.9. 
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Using these linear wave theory estimates as the first guess, the full order optimiza- 
tion, Eq. 6.7, is computed to determine the best full order fit to a global segment of 
the record. The parameters of the potential function computed by this optimization 
are then used as the initial estimate for the final optimization in each defined small 


window in time. 


Phase Shift The phase parameters, a; and a2 are adjusted to accommodate the 


change in the time reference frame from the global to the local solution: 
On = —w,At+ Qn,global (6.23) 


where At is the difference between the time in the two reference frames. 


Nonlinear Optimization in Windows 


The established global solution is used as the first guess for the parameters in 
each window, and these parameters are refined in the same manner as the previous 
chapters by solving Eq. 6.7 with any standard nonlinear optimization routine. 

The resulting solution is then checked to see if is clearly spurious. Spurious solu- 


tions can be identified by the same criteria used in the previous chapters: 
e Very large or systematically variable errors. 
e First order amplitude smaller than one of the higher order amplitudes. 
e Unrealistically large or small frequency or wave number. 
e Large discontinuity between the windows in the predicted kinematics. 


If no solution or a spurious solution is found, the parameters of the solution are 
revised for another attempt. These nevisllons include increasing the width of the 
window, and decreasing the order of the solution. If neither of these adjustments 
result in an acceptable solution, the window is skipped, and future analysis must be 
interpolated through that point. 

As with the analysis of a point pressure measurement, these adjustments are most 


likely to be needed in the long, flat troughs of shallow water waves, or near zero 
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crossings in the record. The additional data provided by the array of measurements 
provides for a more robust optimization than with a single subsurface pressure mea- 
surement. As well as providing additional data, the spatial array provides information 
about the evolution of the wave field in space, helping to define the local wave number 
more clearly. This is in contrast to the analysis of a point measurement, where the 
spatial evolution of the wave field must be determined entirely from the measured 
temporal evolution, coupled with the governing equations. As a result, adjusting the 


solution parameters is necessary less often than with a single measurement. 


6.4 Theoretical Records 


6.4.1 Single Wave Records 


The following results are from “measured” data generated by Fourier steady wave 
theory (Sobey 1989), providing a near-exact solution for steady waves that provides 
the complete kinematics. The data used are a simulation of data that might be col- 
lected by a DWG-1 pressure array (Howell 1992). The DWG-1 is a reliable, easy to 
deploy pressure array. The unit is capable of including up to six independent pressure 
transducers, but is frequently used with three transducers, arranged in an 1.6 m equi- 
lateral triangle (fig. 6.1). Three measurements were chosen for the theoretical data, as 
three is the minimum number required to provide directional information. Additional 
measurements would provide overspecification, and can easily be accommodated in 


the formulation. 


Shallow Water Figure 6.2 is a steady shallow wave generated by 18th order Fourier 
theory with the following parameters: wave height = 3m, period = 10s, water depth 
= 5m, direction of travel 10 degrees from the x-axis, and an opposing Eulerian current 
of 2m/s. This is a near limit wave in this depth water, with this opposing current. 
The pressure array is located at the bottom. The parameters of the LFI solution are: 
window width = ls (7 = 0.1T,), fifth order (J = 5), with 6 samples on the water 


surface, and 6 samples on each of the pressure records (M = J = 6), resulting in 30 
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wave ray 


Pressure Gauges 


Figure 6.1: Layout of DWG-1 pressure array 


equations in 15 unknowns. 

The solid lines are the dynamic pressure, water surface, and kinematics as com- 
puted from the Fourier solution. The circles are the predictions of the LFI method, 
with each circle being located at the center of a separate window. The kinematics 
are measured at a depth of 1.5m below the surface, just under the trough. The LFI 
method has captured the location of the water surface and the kinematics essentially 
exactly, including the pronounced sharp crest. As with the single pressure gauge 
(Sec. 3.4), a fairly high order solution (J = 5) was necessary to capture the sharp 
crest of this shallow water wave. The additional data provided by the three gauges 
allowed for a narrower window than with the single measurement, even at this high 


order. 


Deep Water Figure 6.3 is a steady deep water wave generated by 10th order Fourier 
theory with the following parameters: wave height = 10m, period = 10s, water depth 
= 100m, direction of travel 10 degrees from the x-axis, and zero Eulerian current. The 
pressure array is located 10m below the surface. The parameters of the LFI solution 
are: window width = 1s (7 = 0.1T.), second order(J = 2), with 3 samples on the 


water surface, and 3 samples on each of the pressure records (M = J = 3), resulting 


Kinematics at the center of a pressure array 


— Fourier 
Oo LFI 


tw. 


Figure 6.2: LFI predictions (J = 5) and exact kinematics at the center of the array 


at z = —1.5 m for a steady shallow water wave. 
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Kinematics at the center of a pressure array 


Figure 6.3: LFI predictions (J = 2) and exact kinematics at the center of the array 


at z = —5 m for a steady deep water wave. 
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in 12 equations in 9 unknowns. 

The LFI method has captured the location of the water surface and the kinematics 
essentially exactly, at only second order. Higher order may be necessary to define the 
kinematics of a more irregular record, and can be included with little computational 


penalty. 


6.4.2 Records of Two Intersecting Waves 


Theoretical pressure records were generated by the same fourth order stokes-type 
intersecting wave theory used in the previous chapter (Ohyama, Jeng, and Hsu 1995a). 
This method provides the complete kinematics, is applicable in deep water, and as- 


sumes a zero Eulerian current (see appendix A). 


Standing Wave Figure 6.4 shows the results of the method for a standing wave. 
The parameters of the Ohyama solution are: period = 10s, first order amplitudes = 
3m (resulting in a total wave height of slightly over 6m), propagation directions = 15 
and 195 degrees from the x axis, in deep water. The pressure array is located 10m 
below the mean water level. The LFI method applied to third order (J = 3) finds 
the the water surface and the kinematics at the elevation of 3m below the surface, 
just below the trough of the wave. While these results show the complete wave, it is 
important to keep in mind that as with all the previous results, each of the indicated 
points is in the center of a separate window, and was computed independently of the 
other windows. In this case, the standard window width was 1s (7 = 0.1T.), with 
a third order potential function (J = 3), four water surface nodes (M = 4) and six 
samples on the pressure records (J = 6) distributed equally in time in each window, 


resulting in 26 equations in 24 unknowns. 


Short Crested Wave Figure 6.6 shows the results of the method for the short 
crested wave shown in Figure 6.5 The parameters of the wave are: period = 10s, 
first order amplitudes = 3m (resulting in a wave height of just over 6m), propagation 


directions: 30 and 150 degrees from the x axis, in deep water. The pressure array 
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Kinematics at the center of a pressure array 


Figure 6.4: LFI predictions (J = 3) and exact kinematics at the center of the array 


at z = —3 m for a standing wave. 
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300 


Figure 6.5: Water surface of short crested wave at t=0 


is located at 10m below the MWL. The LFI method applied to third order (J = 3) 
again finds the kinematics for this wave almost exactly. In this case, the standard 
window width was also ls (7) = 0.17.), with four water surface nodes (M = 4) and 
six samples on the pressure records (J = 6) distributed uniformly in time in each 
window, resulting in 26 equations in 24 unknowns. 

Figure 6.7 is the same short crested wave, but computed with the single wave 
method of the LFI solution. In this case, the LFI solution still matches the measured 
pressure record very well, as is virtually always the case, as those measurements are 
part of the optimization. The predicted water surface is also fairly accurate, but 
with the predicted crest slightly underestimated. The vertical velocity is also fairly 
accurate. The LFI predictions for horizontal velocities, on the other hand, are very 
different from the Ohyama solution. The resulting discontinuities in the predictions 
for the horizontal velocities make it clear that something important is missing from 
the solution. The single wave method is simply not able to capture a short crested 


sea made up of two distinct directional components. 


Kinematics at the center of a pressure array 
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Figure 6.6: LFI predictions (J = 3) and exact kinematics at the center of the array 


at z = —3 m for a short crested wave. 
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Kinematics at the center of a pressure array 


la 


Figure 6.7: LFI predictions (J = 3) and exact kinematics at the center of the array 
at z = —3 m for a short crested wave computed with the single wave method. 
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Choice of Order 


The choice of order for the analysis of arrays of pressure measurements is essen- 
tially the same as for a single pressure measurement. See chapter 3 for a detailed 
discussion. In general, lower orders are necessary for deep water than shallow water, 
with higher than third order unlikely to be necessary in deep water. Shallow water 


may require up to fifth or sixth order in order to capture near limit waves. 


Choice of Window Width 


The choice of window width is also very similar for the analysis of pressure arrays 
as it is for the analysis of a single pressure measurement, or arrays of water surface 
measurements. Window widths of a minimum of about 1/10 the zero crossing period 
are required, with the occasional need to widen the windows near zero crossings. The 
additional data provided by the array of pressure measurements allowed a solution to 
a very steep shallow water wave at high order with a window width of 0.17,, where 
a similar solution required the window width to be doubled to 0.27, when there was 


only a single point pressure measurement available (Section 3.4). 


Number of Nodes on the Water Surface and Pressure Records 


The criteria developed in the previous chapters for the number of water surface 
nodes are applicable here. J+1 points are required to specify uniquely a Fourier series 
of order J. In order to keep the order of the water surface consistent with the order 
of the potential function, at least M = J+1 water surface nodes should be used. The 
same number of points on the measured pressure records is usually adequate. For the 
standing and short-crested waves, additional equations were needed to specify the 
solution, so J = J + 3 points on the pressure records were used. The results are not 
highly sensitive to this parameter, provided sufficient samples are used to define the 
local curvature and specify the system of equations. More points on the measured 


pressure records will assist in accommodating noise in the record. 
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6.5 Field Measurements 


The following results are from data provided by Gary L. Howell of the US Army 
Corps of Engineers Waterways Experiment Station, Coastal Engineering Research 
Center (CERC). The data were collected as part of CERC’s Coastal Inlets Research 
Program. They were collected near the Ponce de Leon inlet, Florida, by a DWG-1 
with three absolute pressure gauges, data sampled at 5 Hz. The DWG-1 was resting 
on the bottom, so that the pressure sensor diaphragms were positioned 0.21m from 
the bed. 


6.5.1 Field Results 


Figure 6.8 is a segment of a record collected on August 20, 1996. The statistics of 
the record, based on about one hour of record are: mean wave height = 1.05m, peak 
period = 5.6s, peak direction = 83° from true north, in a depth of approximately 
7.9m. 

The top plot is the dynamic pressure at the center of the array. The solid line 
is the approximate measured pressure, computed by linear interpolation from the 
three measured points in the array. The circles are the values predicted by the 
LFI method. The other four plots are the predicted water surface, and the three 
orthogonal velocities, as predicted by the LFI computation. The single wave method 
was used, with the following parameters: third order (J = 3), four water surface 
nodes and four samples on the pressure records (J = M = 4), and a window width of 
1.42s (0.37-), resulting in 20 equations in 11 unknowns. There were no data available 


about the Eulerian current, so it is taken as zero. 


Window Width In the previous section, results on theoretically generated records . 
indicated the successful solutions were possible with quite narrow windows, gener- 
ally one tenth of the zero crossing period. When working with this field record, the 
optimization converged with a window width that narrow, but it resulted in wild fluc- 
tuations in the predicted propagation direction from window to window. Figure 6.9 is 


the same segment of the record as in figure 6.8, computed with the same parameters, 
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Predicted kinematics at the center of the array 


Figure 6.8: LFI predictions (J = 3) of kinematics at the center of the array at the 


water surface for a field record. 
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except that the window width used was much narrower, 7 = 0.1T7-.). There are large 
discontinuities in the horizontal velocity in the z direction (third plot), near both 
crests. Figure 6.10 shows the parameters of the potential function from this solution. 
6, where 0 = tan7!(k,/kz), is the local direction of propagation of the wave. There is 
a sudden shift in direction of 80° at the first crest, and over 100° at the second crest. 
This would result in unrealistic predictions of huge accelerations. Figure 6.11 shows 
the parameters of the solution computed with 7) = 0.37. With the wider window, 
the propagation direction varies smoothly around the peak direction of the record, 
Q = 83°. 


Array Size The need for a fairly wide window for the field solution is likely due 
to the small size of the pressure gauge array. Compact size is one of the major 
advantages of the DWG-1 array, as it is a single unit that is easily deployed, but 
it makes it susceptible to difficulties with local analysis. The propagation direction 
within each window is determined by the spatial gradients of the dynamic pressure, 
as estimated by the measurements. In segments of the record where the gradients 
are small, the predicted direction of propagation can fluctuate a great deal with only 
small change in any of the measured values. These changes could be the result of 
measurement error, or be small fluctuations in the actual dynamic pressure. In either 
case, the large resulting change in predicted propagation direction can be a source of 
error. 

The linear dispersion relation predicts that the wavelength corresponding to the 
peak period of this record in this depth of water would be over 42 m. The DWG-1 
array is only 1.6 mon a side, which is about 0.04 of the approximate wavelength of the 
peak period waves. As a result, near the crest of the wave, where the water surface 
is close to horizontal, all three gauges are near this point, and the local estimate 
of the gradient of the pressure is very sensitive to small fluctuations. This effect 
is mitigated by using a larger window in time, so that a part of the wave with a 
higher gradient is part of the window, allowing the propagation direction to be better 
defined, and smoothing out any measurement errors. As minimum window size has 


been determined to be about 1/10 of the zero crossing frequency, an array size of 
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Predicted kinematics at the center of the array 
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Figure 6.9: LFI predictions (J = 3) of kinematics at the center of the array at the 


water surface for a field record, using a narrow window. 
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Solution Parameters 
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Figure 6.10: Parameters of the LFI solution for a field record, using a narrow window 
(To = UEIIE 
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Solution Parameters 


nat 


Figure 6.11: Parameters of the LFI solution for a field record, using a wide window 
(To = Orsi): 
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about 1/10 of a zero crossing wavelength would probably result in more stable results 


with the smaller windows. 


6.6 Discussion 


The LFI method can be applied to interpretation of records from arrays of pressure 
measurements. It is effective for records from primarily unimodal seas, as well as the 
bimodal seas that are generated near reflecting surfaces. The result of the analysis is 
a complete description of the kinematics of the waves throughout the water depth, in 
the vicinity of the array. The method is local in that a separate potential function is 
used to describe each small segment of the record, using information from only that 
segment. This approach is rational in that it seeks to satisfy the governing equations 
of gravity waves, including the full nonlinear free surface boundary conditions. 

Using a potential function representing a single progressive wave, the method 
has been shown to be effective on theoretically generated records of steady waves in 
both deep and shallow water. This formulation is appropriate in locations far from 
reflecting surfaces, where the sea state is expected to be unimodal. When the method 
is applied with a potential function representing two intersecting waves, it is effective 
' in accurately capturing the kinematics of standing and short crested waves, as result 
from the reflection of steady waves from a reflecting surface, such as a sea wall. This 
two wave formulation is appropriate near such reflecting surfaces, where the sea state 
is likely to be bimodal. 

A number of parameters must be specified to establish a solution, including the 
order, window width, number of water surface nodes defined, and number of samples 
on the measured records. The examples in this chapter indicate that the criteria 
used to define these parameters are similar to those developed in the analysis of point 
pressure records and arrays of water surface measurements. Fairly low order solutions 
(J = 2 or 3) are effective in deep water, while shallow water requires higher order 
solutions (J = 5 or 6). 

Window widths of 1/10 the zero crossing period of the record are effective on 


theoretical records. With field records, however, when the array used is small (less 


152 


than 1/10 a typical wavelength), the predictions of the propagation direction can 
be very sensitive to subtle fluctuations in the measure record. In order to find a 
stable solution, the window width must be increased to about 1/3 of the zero crossing 
frequency. This includes more of the record in the window, allowing the propagation 
direction to be more clearly defined. 

The examples in this Chapter were all computed using MATLAB running under 
the Linux operating system on an Intel 90MHz Pentium processor based PC. The 
shallow water steady wave (Fig 6.2) took the longest to compute, about 64min., and 
471x10° floating point operations (flops). The steady deep water wave (Fig. 6.3) took 
the least time, 12min. and 22x10° flops, the standing wave (Fig. 6.4) took 16min. 
and 131x10® flops, the short crested wave (Fig. 6.6) took 13min. and 133x10° flops, 
and the field record (Fig. 6.8) took 17min. and 46 x10® flops. 
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Chapter 7 
Conclusions 


This dissertation examines the problem of determining the kinematics of three di- 
mensional irregular seas from a variety of wave measurements. The presented methods 


and results lead to the following conclusions: 


e Currently used global methods are inadequate for determining accurately the 


kinematics of irregular waves. 


— Linear superposition does not satisfy the nonlinear free surface boundary 
conditions, resulting in large errors in the predicted kinematics near the 


water surface. 


— Directional spectral methods generally disregard the phase information, 
making it impossible to determine the detailed kinematics of directional 


seas. 


— Global methods of high enough order to accurately solve the free surface 
boundary conditions in three dimensions are prohibitively computationally 


intensive. 


e The Local Fourier Method for Irregular waves (LFI) can accurately describe the 
kinematics of two and three dimensional irregular seas by satisfying the full free 


surface boundary conditions. 
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— The local nature of the LFI method allows it to satisfy the nonlinear free 
surface boundary conditions at low order by finding separate solutions in 


each of a sequence of small windows in time. 


— The LFI method provides an accurate two dimensional description of the 
kinematics of irregular seas measured by a single water surface gauge or 


subsurface pressure gauge. 


— The LFI method provides an accurate three dimensional description of the 
kinematics of irregular seas measured by arrays of water surface gauges or 
arrays of subsurface pressure gauges, including the bimodal seas resulting 


from reflection from a vertical surface, such as a sea-wall or breakwater. 


— When working with subsurface pressure records, the LFI method is limited 
in its capability to capture the high frequency components of the sea state 


that decay rapidly with depth. 


— In order to accurately capture the detail of the measured records, the LFI 


method requires very accurate data, with high sampling rates. 


— The LFI method can be adapted to virtually any arrangement of wave 


measuring instruments. 


— The LFI method required substantial, but not prohibitive, computational 


resources. 


7.1 Future Work 


The results presented indicate the promise of the LFI method, and local methods 
in general. There is still much work to be done before the method could be considered 


” code. 


usable as a “black box’ 

Much of this work revolves around determining the numerical details of the nonlin- 
ear optimization. While the presented results provide some guidelines to determining 
appropriate parameters, considerable judgment and experimentation with any given 
record is required. In the future, it may be possible to establish criteria for defining 


the following numerical solution parameters: 
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e Window width 
e Order of the potential function 
e Number of water surface nodes 


Number of samples of the measured record 


e Number of iterations allowed in the optimization 
Some of the criteria that might be used to help define these values are: 
e Non-dimensional depth 
e Ursell number 
e Local curvature of record 


In principle, the LFI method may be extended to additional arrangements of in- 
struments. While it is straightforward, in theory, to adapt the method to virtually 
any arrangement of instruments, each different arrangement is likely to present a new 
set of numerical solution difficulties. One of the sources of this is the use of non- 
_ linear optimization. Different optimization problems tend to be unique, and require 
substantial experience in order to determine the appropriate approach. 

The experience gained in applying the LFI method to a variety of arrangements 
of instruments helps demonstrate the need for appropriate data. Whether this par- 
ticular method of analysis is employed, or any other nonlinear technique, the need for 
accurate data is paramount. The design of a field data collection program inherently 
includes assumptions as to how the resulting data will be analyzed. Currently used 
data collection programs are often designed with the assumption that the data will be 
analyzed with the common methods of linear spectral analysis. As a result the data 
may not be appropriate for the nonlinear analysis needed to determine the detailed 
kinematics of a measured wave field. 

Some of the considerations that should be included in the design of a field mea- 


surement program are: 
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e Array size appropriate for the analysis method chosen and wave field expected. 
e Measurement or prediction of the local current field. 


e Sampling rate that captures the detail of the record. 


Reduction of sensor noise and increased instrument accuracy. 


As the methods for measuring waves and interpreting those measurements become 
more accurate, the understanding of the complex processes in the ocean and coasts 
that are molded by waves will increase as well. This understanding should lead to safer 
coastal and offshore structures, and well as reducing the impact that these structures 


have on the sensitive coastal environment. 
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Appendix A 


Intersecting wave theory 


A.1 Introduction 


The theory of multiple wave interactions presented by Ohyama, Jeng, and Hsu 
(1995a) has been used for two purposes in this work. The most important function 
was to suggest a form for a general three dimensional potential function representing 
interacting waves. The form chosen is presented in Chapter 4. It also provided a 
theoretical method for generating a complete, consistent set of records to use for the 
development and testing of the LFI method. This appendix presents a review and 


detailed critique of the intersecting wave theory. 


A.2 Theoretical Background 


Ohyama et al.’s method is an extension of well accepted and verified methods. 
It is essentially a Stokes type expansion in the wave steepness. As such, it is ex- 
pected have greatest applicability in deep water. The flow is taken to be irrotational 
and incompressible, and thus the governing equations are the same as those used in 


chapter 4 and for most finite depth irrotational wave theory: 
e Mass conservation (Laplace equation) (Eq. 4.2) 


e Bottom boundary condition (Eq. 4.3) 
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e Dynamic free surface boundary condition (Eq. 4.4) 
e Modified kinematic free surface boundary condition (Eq. 4.6) 


There is no accommodation for an Eulerian current, although it might easily be 
included in future work. The non-dimensional potential function, water surface, and 
frequencies are expanded in a power series in the small parameter, (€ = ka), where k 
is a typical wave number and a is a typical amplitude of the first order component of 


the water surface. 


b=) + 2d + Hd) + AG + O(€) (A.1) 
p= GB Reg? + eg” + Gg £0) (A.2) 
w= eu”) + ew?) + eu” + Au) + O(e?) (A.3) 


Where the dé”, n'” are the potential function and water surface at order n, and ual”) is 


the frequency of the zth wave at order n. The steepness of all the intersecting waves 
is taken to be of the same order. Modern Stokes theories have found it necessary 
to use an expansion parameter based on the physical wave height, rather than the 
amplitude of the first order component (Fenton 1990). This is not feasible with 
intersecting waves, as there is no clearly defined wave height. While not ideal, Stokes 
methods based on this first order expansion parameter can be effective (Skjelbreia 
and Hendrickson 1962), as long as they are algebraically correct and not pushed to 
near limit waves or shallow water. 

Care must be taken when computing a given order solution when multiple parts 
of the solution are all expanded in a perturbation series. In this case, the potential 
function, water surface, and frequencies are expanded separately, but the potential 
function and water surface are dependent on the full order frequency. Ideally, a 
given order component will depend only upon parameters of the same order. In this 
case, however, the full order frequency must first be computed, and then used as 
the frequency at all orders. If the order of the frequency used was matched to the 


order of 7 or ¢ being computed, the resulting solutions would be out of phase, with, 
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for example, 6") having a different frequency to 6). It could be argued that this 
approach would be accurate to second order, but the components would be shifted 
out of phase as time progressed. It is more accurate to compute the results using the 
full order frequencies at all orders of ¢ and 7. 

The first order solution for N intersecting waves is the familiar superposition of 


linear waves: 


N 
(@) 2 a; cosh k(h + 2) a (if ree 4 
ro) De i a eee sin (k,z + kyy — wt + a) (A.4) 
7) = yy a; cos (ky + kyy — wt + a) (A.5) 


i=1 


w'') = ./gk; tanh kh; (A.6) 


The higher order solutions are computed by applying the Laplace equation, the 
bottom boundary condition, and the free surface boundary conditions, expanded in 
Taylor series about the mean water level. The algebra and the resulting solutions 
are long and extremely tedious, are presented in Ohyama et al. (1995a), and they 
will not be repeated here. It should be noted that the frequency was only solved to 
third order, as the fourth order frequency would require a solution to the fifth order 


potential function, which was not presented. 


A.3 Analytical Verification 


Ohyama et al. provided a number of comparisons with other wave theories. When 
the method is applied for a single wave component, it represents a steady wave, and as 
such should be expected to provide the same solution as conventional Stokes theory. 
The authors presented a comparison with the Fenton (1985a) solution. The Fenton 
solution is based on an expansion parameter based on the wave height, rather than 
the first order amplitude: € = kH/2. By solving for € in terms of the «, The two 


solutions were found to coincide exactly. 
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When applied for two wave components of the same wavelength, same amplitude, 
and opposing directions, the solution is a standing wave. The resulting solution 
coincides almost exactly with a standing wave solution given by Chen (1988), except 
for one coefficient. The authors suggest that the slight discrepancy is the result of a 
typographical error in Chen’s paper. 

Ohyama et al. also indicated that they compared their solution with those for 
short-crested seas (Hsu 1990; Hsu and Chen 1992; Fenton 1985c) and once again 


found essentially perfect agreement. 


A.4 Numerical Verification 


The above theoretical verifications indicate that the method is correct when re- 
duced to two particular simplifications. It remains to verify that the method is 
accurate when applied to a more complex situation, with a number of intersecting 
waves of different heights, directions and wavelengths. 

The actual computation of the general solution is quite complicated and contains 
a very large number of coefficients that must be computed. Despite this complexity, 
once the code is written and debugged, it is a trivial matter for modern computers 
- to calculate the full solution. It is not, however, a trivial matter to write the code 
without errors or, indeed, to find the likely typographical errors in the published 
paper. A number of these typographical errors have been presented in the published 
erratum (Ohyama et al. 1995b). 

An additional error has been found and confirmed by personal communication 


with the authors. pp. 55 Eq. 47 in Ohyama, Jeng, and Hsu (1995a) should read: 


¢ = {£?Bo20 = Baro} cos 2kx + €4Bs49 cos 4kr + {2 2s "Baus } cos kx cos wt 
+ ror + <“Bixy} cos 2kz cos 2wt + €°Ba13 cos kx cos 3wt 
SES hee cos 3kz coswt + €? Boga cos 3kz cos 3wt + a" Bana cos 2kz cos 4wt 


+€4B442 cos 4k cos Qwt + €4B444 cos 4k cos 4wt + Ofe*} 


As a result of the complexity of the computational code, it is advisable to use numer- 
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ical methods to verify the validity of the solution and resulting computational code. 
The following results were produced by a Fortran code derived from the code written 
and generously provided by Takumi Ohyama that was used to compute the water 
surface plots given in (Ohyama et al. 1995a). The code was expanded to compute 


the full kinematics. 


A.4.1 Richardson Extrapolation 


Fenton (1985a) presented a variation of Richardson extrapolation that can be used 
for numerical verification of wave theory code. Richardson extrapolation (also known 
as extrapolation to the limit) is a method traditionally used to increase the accuracy 
of finite difference computations (Salvadori and Baron 1961; Roach 1976). In Fenton’s 
adaptation, the method is used to simultaneously check all the terms in the solution, 
and to provide an estimate of the order of the accuracy of the result. 

The form of Ohyama et al.’s solution solves the field equation (Laplace) and the 
bottom boundary condition exactly. It is expected to satisfy the free surface boundary 
conditions to the order it is applied. The following method calculates the order of 
accuracy of the free surface boundary conditions. Classical Richardson extrapolation 
is used to evaluate the error of a solution at a given point in time and space. By 
taking advantage of the periodicity of wave motion, Fenton’s variation can be used 


to evaluate the solution throughout a complete wave. 


Single Wave 


The solution for a single wave is periodic and symmetric about the crest, and thus 
the error in each of the boundary conditions can be expanded in a single sided Fourier 


series: 


ges ge” (A.7) 
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[ Fourier Component, (| 0] 1)2)3]4]5] 6] 
[_ Order of KFSBC error [2.0 [2.0/2.0] 2.02.01 2.0] 2.0] 
[_Order of DFSBC error [3.0] 2.2 [2.0 [2.02.0] 2.0] 200 | 


Table A.1: Order of errors for a single wave computed to order 1 


[Fourier Component, (K)] 0] 1]2]3]4][5] 6 
[_ Order of KFSBC error | 3.0/3.0 30 [30/30] 3.0/3.0] 
[Order of DFSBC enor [2.9/3.0] 3.0/3.0] 3.0/3.0] 3.0] 


Table A.2: Order of errors for a single wave computed to order 2 


The magnitude of this error, and thus each of the coefficients can be expanded in a 


first order Taylor series. 


Cy(Ex) = Bub. + O (&) (A.8) 


These coefficients are a function of the small parameter, € and thus are a function of 


€” at nth order, &, =e". 
C.(€) = Bre +O (e*) (A.9) 


When the discreet transforms of the errors are computed for two different values of €, 
an approximation for the n, can be computed, giving an approximation to the order 
of the errors of each harmonic in each boundary condition. 


iy, = GaGa + O(€1, €2) (A.10) 


Table A.1 presents the results for a single wave, computed to order | in deep water 
(kh = 100, & = 0.01, €2 = 0.02). These results show that the error is of at least 
second order in both of the boundary conditions, indicating that the computation is 
accurate to first order. 

Tables A.2 through A.4 present the results for the same wave, computed to sec- 


ond through fourth order. These results indicate that the error in both boundary 
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[Fourier Component, ()] 0] 1]2]3]4]5] 6 
[_ Order of KFSBC error [4.0 [4.0] 4.0 [4.0 [4.0 [4.0/4.0] 
[Order of DFSBC error | 4.0/3.9] 4.0] 41 [4.0] 40] 40] 


Table A.3: Order of errors for a single wave computed to order 3 


[Fourier Component, (K)] 0] 1]2[3]4]5] 6] 
[Order of KFSBC eor [5.0 [5.0 [5.0] 5.0] 5.0 |5.0 | 5.0 | 
[_ Order of DFSBC error [48 [5.0/5.0] 5.0 [5.0 [5.050] 


Table A.4: Order of errors for a single wave computed to order 4 


conditions is of at least one order higher than the order used to compute the results. 
These results serve to confirm the efficacy of the Richardson extrapolation method, 


as well as the accuracy of the wave theory up to fourth order. 


Standing Wave 


Ohyama et al.’s method can also be used to compute a standing wave by defining 
two waves of the same wavelength and amplitude, and opposing directions. In this 
case, the wave is not steady, and the errors in the free surface boundary conditions 
must be considered for a full period in time and full wavelength in space. The ex- 
trapolation to the limit method can be extended to this case by expanding the free 


surface errors in a two dimensional Fourier series (Newland 1993). 


=) Ss Cp (jee e/E+H/7) (A.11) 


k=0 I=0 

An estimate of the order of the errors can then be computed in a similar manner as Eq. 
A.10. The results of this computation for a standing wave in deep water (kh = 100, 
€; = 0.01, €2 = 0.02) are given in Tables A.5 and A.6. There are no values given 
for the zeroth order in time of the kinematic boundary condition because C(€),,o for 
both values of € are zero to the precision of the calculation, so the values computed 


are spurious. These results indicate that Ohyama et al.’s method is also accurate to 
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[Fourier Component, [0] 1]2[3]4)]5] 6] 
(ea 
[1 [50 [50 [50 [50 [5050/50 
en dd Ne 
[Es 5050/8050 [5.05.0 | 5 | 
a 
[ES 49 [51 [as [50 [51 50) 5.0 
[Ee [5050 [51 [5.0 [5.0 [50] 5.3 | 


Table A.5: Order of KFSBC errors for a standing wave computed to order 4 


[Fourier Component. (0 ]1]2[/3]4]5] 6] 
ee EES 
a4 [50/4150] 49]50 [50 | 


[61 [52 [6.5 [5.0] 5.8 [5.5 [5.3 | 
[5 [495.0 [5.2] 5152] 8.049 | 
[E853 [5.0] 5.95.0] 47 [8.0 | 53 | 


Table A.6: Order of DFSBC errors for a standing wave computed to order 4 


approximately fifth order for a standing wave. 


A.4.2 More Complex Seas 


Fenton’s method relies on the periodicity of the solution to be tested. If the 
solution is not periodic in space and time, then the errors in the boundary conditions 
will not be periodic, and the Fourier expansion cannot be strictly applied. Fourier 
coefficients could still be computed by assuming periodicity, but the discontinuity 
generated by the assumption of periodicity would generate spurious results. These 
could be minimized by taking a large number of points over a large range in time and 
space, but this would become very computationally intensive. 


The accuracy of any approximate solution can be measured by examining the 
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oder (a) 
342 x 107 


Table A.7: RMS errors in free surface boundary conditions for a short-crested sea 


errors in the governing equations resulting from the solution. Fenton’s method is a 
particularly powerful method for examining these errors, in that it considers an entire 
wavelength and period at once and gives an estimation of the order of accuracy of the 
solution. A simpler approach can also be useful and be universally applied. In the 
case of Ohyama et al.’s method, the field equation and bottom boundary condition 
are exactly satisfied at all orders of approximation. If the solution is accurate, the 
errors in the free surface boundary conditions will decrease with each increase in the 
order of approximation. While this method does not guarantee that the solution is 
correct, most errors in the algebra or the computational code will result in a decrease 
in the accuracy of the result at the order in which the error occurs. 

Table A.7 contains the errors in the free surface boundary conditions at orders 
one through four for a short crested sea. This sea is created by two waves of the same 
height and wavelength (« = 0.1, kh = 100) intersecting at an angle of 30 degrees. 
The values given are the root mean square of the errors computed by averaging over 
a grid of points one wavelength in both horizontal directions, and over one period. In 
this case the errors in the free surface boundary conditions decrease with increasing 
order of solution. It should be noted, however, that the errors do not decrease as 
much as e”. This might be partially due to the fact that the frequency is expressed to 
only second order. Including the fourth order frequency would probably reduce the 
error at fourth order. This behavior warrants closer examination, particularly if the 


method were to be expanded to higher order. 
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A.5 Depths of Validity 


As the method is very similar to Stokes theory, it is expected to be valid in a similar 
range of depths. Stokes theory is valid if both ¢ and the parameter, ¢/(kh)*, similar 
to the Ursell number, are small (Fenton 1985a). This indicates that the method has 


optimal applicability in deep water, and should not be applied in shallow water. 
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